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Abstract 

This  investigation  develops  an  interactive  computer 
aided  method  for  designing  lowpass  or  highpass  two- 
dimensional  finite  impulse  response  digital  filters.  Filters 
designed  using  this  method  will  have  linear  phase  and  an 
equiripple  error  in  the  transmission  and  attenuation  bands. 

The  essence  of  the  method  is  that  it  transforms  an 
optimal  one-dimensional  digital  filter  into  a close  approxi- 
mation of  an  optimal  two-dimensional  digital  filter.  The 
amplitude  characteristic  of  the  one-dimensional  filter 
is  preserved  in  the  sense  that  each  point  of  the  one- 
dimensional frequency  response  is  mapped  to  a contour  in 
the  two-dimensional  plane.  This  transformation  was  first 
proposed  by  James  H.  McClellan,  and  is  now  called  McClellan 
Transformation . 

By  controlling  the  mapping  of  a specified  one- 
dimensional frequency  to  a desired  contour  shape  in  the 
plane,  two-dimensional  filters  of  fairly  arbitrary  specifi- 
cations can  be  designed;  that  is,  their  frequency  response 
can  be  determined,  and  the  associated  two-dimensional 
impulse  response  coefficients  calculated. 
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INVESTIGATION  OF  OPTIMAL 
LINEAR  SHIFT- INVARIANT 
TWO-DIMENSIONAL  DIGITAL  FILTERS 

I.  Introduction 

This  investigation  concerns  an  interactive  computer 
aided  method  for  designing  two-dimensional  finite  impulse 
response  (FIR)  digital  filters.  A spectral  transformation 
technique,  suggested  by  McClellan  (Ref  1),  is  extended  to 
permit  the  design  of  two-dimensional  filters  within  a 
special  class.  Through  this  transformation,  an  optimal  one- 
dimensional FIR  filter  is  transformed  into  an  equiripple, 
though  sub-optimal,  two-dimensional  digital  filter.  The 
amplitude  characteristic  of  the  one-dimensional  filter  is 
preserved  in  the  sense  that  each  point  of  the  one-dimensional 
frequency  response  is  mapped  to  a contour  in  the  two- 
dimensional  frequency  plane. 

Digital  Filters 

Digital  systems  deal  with  signals  that  are  discrete  in 
both  amplitude  and  time.  Digital  filtering,  sometimes  called 
digital  signal  processing,  is  involved  with  operations  on 
such  signals,  using  either  general  purpose  computers  or 
special  purpose  hardware.  In  the  early  years  of  digital 
filtering,  the  required  processing  could  not  always  be  done 
in  real  time.  However,  advances  such  as  the  Fast  Fourier 
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Transform  algorithm  now  permit  processing  times  several 
orders  of  magnitude  faster  (Ref  2:3).  Although  digital 
filters  are  often  used  as  mere  simulations  of  analog  filters, 
their  most  significant  value  is  when  they  are  used  to  create 
filters  with  characteristics  unobtainable  using  analog 
methods . 

To  permit  the  use  of  well-known  linear  system  theory, 
only  linear  shift-invariant  (LSI)  digital  filters  will  be 
considered  in  this  investigation.  LSI  filters  can  be 
divided  into  two  classes  depending  on  the  duration  of  the 
filter's  impulse  response;  infinite  impulse  response  (HR) 
filters  (or  recursive  filters),  and  finite  impulse  response 
(FIR)  filters  (or  nonrecursive  filters).  IIR  filters  were 
originally  more  popular  than  FIR  filters,  because  the 
traditional  approach  to  digital  filter  design  involves  the 
transformation  of  an  analog  IIR  filter  into  a digital  IIR 
filter  with  some  prescribed  specification  (Ref  2:197). 
However,  FIR  filters  have  two  attractive  properties.  First, 
there  is  the  possibility  of  designing  exact  linear  phase, 
required  in  many  applications  (Ref  7*76).  Second,  the 
FIR  filter  is  never  unstable.  These  qualitites  are  often 
important  in  digital  signal  processing  applications.  For 
many  years,  efficient  methods  for  designing  optimal  1-D  FIR 
filters  were  not  available;  however,  in  1972,  Parks  and 
McClellan  published  a very  efficient  algorithm  for  the 
design  of  optimal  1-D  FIR  filters  (Ref  3»  Ref  4). 

With  this  problem  solved,  attention  turned  to  the  design 
of  filters  of  higher  dimensions,  especially  two-dimensional 
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(2-D)  digital  filters.  2-D  filters  find  applications  in 
areas  where  signal  specification  requires  spatial  coordi- 
nates, such  as  in  image  processing  (Ref  5)  or  seismic 
analysis  (Ref  6).  A digital  filter  must  be  designed  (that 
is,  the  frequency  response  must  be  specified  and  the  impulse 
response  coefficients  must  be  determined) , and  then  the 
filter  is  implemented.  To  permit  the  proper  amount  of  detail 
in  this  investigation,  only  the  design  of  the  digital  filter 
will  be  considered.  Recognizing  that  implementation  is 
certainly  just  as  important  as  design,  a short  listing  of 
implementation  references  has  been  included  as  a supplemen- 
tary bibliography. 

The  Design  Problem  of  2-D  Filters 

A design  goal  for  a 1-D  filter  might  be  to  minimize 
the  error  over  the  interval  rr , rr^  between  a desired 
frequency  response  D(w)  and  the  obtainable  frequency  re- 
sponse H(cj),  according  to  some  error  criteria.  It  is 
assumed  that  the  frequencies  of  interest  have  been  normal- 
ized to  the  interval  ^-Tr.trJ.  This  is  the  case  if  the 
sampling  time  T is  equal  to  unity.  Alternatively,  one 
maps  the  sampled  frequencies  via  the  complex  exponential 
mapping  z = exp(sT)  and  deals  with  the  set  of  angles,  repeating 
£-Tr,TrJ  on  the  Z-plane  unit  circle  (Ref  2*199-200). 

The  design  goal  in  two -dimens ions  is  analogous,  where 
the  interval  [-tt.ttJ  is  replaced  by  the  region  £-rr,TrJ  x[-Tr,Trl 


3 


in  the  spatial  plane.  ' 2-D  FIR  filters  can  be  designed  by 
multiplying  the  infinite-extent  ideal  frequency  response  by 
an  NxN  sample  point  rectangular  window.  This  amounts  to  a 
truncation  of  the  infinite  duration  impulse  response  sequence. 
However,  as  in  1-D,  the  rectangular  window  produces  a large, 
undesirable  ripple  in  the  frequency  response,  due  to  the 
Gibb's  effect.  By  using  better  windows,  the  ripple  can  be 
reduced  somewhat,  but  the  resulting  filter  is  never  optimal, 
since  the  required  convolution  "smears"  the  frequency 
response  (Ref  7:239-250). 

Rabiner  and  others  have  designed  optimal  2-D  FIR  filters 
by  using  linear  programming.  While  the  resulting  filters  were 
excellent,  the  design  method  was  very  costly  and  inefficient. 
The  highest  order  filter  designed  (9x9  sample  points)  in- 
volved thousands  of  constraint  equations  and  required  more 
them  one  hour  of  computation  time  on  a high  speed  computer 
system  *(Ref  ?«471). 

Effort  has  been  made  to  generalize  McClellan  and  Park's 
1-D  algorithm  to  higher  dimensions  despite  an  assertion  by 
McClellan  that  the  algorithm  could  not  be  extended  directly 
(Ref  1»30).  Both  Kamp  and  Thiran  (Ref  8)  and  Hersey  and 
Mercereau  (Ref  9)  recently  published  algorithms  for  the 
design  of  optimal  2-D  FIR  filters  based  on  variations  of  the 
Parks-McClellan  method.  Despite  use  of  efficient  numerical 
techniques,  no  dramatic  decrease  in  execution  time  was 
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achieved  in  comparison  with  linear  programming.  Addition- 
ally, these  methods  are  restricted  to  low  order  filters 
(Ref  10s 240). 

In  1973 » McClellan  proposed  that  a more  efficient  method 
of -designing  optimal  2-D  filters  would  be  by  a "spectral 
transformation" , a complex  mapping  that  carries  a stable 
rational  transfer  function  into  another  stable  rational 
transfer  function.  The  change  of  variables  McClellan  sug- 
gested, to  be  referred  to  as  the  McClellan  Transformation,  is 

cosu  = F(w^,to2)  = t^  + t2coso>1  + t^cosw2  + t^cosco^  * coso>2  (1) 

Using  this  transformation,  with  the  values  t-^  = -.5, 

^2  = ^3  = ^4  = McClellan  was  able  to  design  circularly 

symmetric  2-D  filters  by  transforming  a suitable  1-D  optimal 
filter  (Ref  1).  The  magnitude  of  each  frequency  in  the 
interval  [j-rr , tt J is  mapped  to  an  4d^-g>2  contour  in  the  square 
region  [-rr,  rrjx  J^-rr , -rrj  . Several  such  contours  are  shown  in 
Fig.  1,  and  the  resulting  2-D  frequency  response  in  Fig.  2. 
(The  design  program  developed  in  this  investigation  was 
used  to  obtain  these  figures.)  Because  there  are  many  2-D 
applications  that  either  require  or  can  tolerate  a 
circularly  symmetric  filter,  McClellan’s  method  has  been 
widely  praised  (Ref  9*1-2). 

Because  these  filters  are  obtained  using  a transform 
method,  they  are  generally  sub-optimal.  However,  they  do 
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Fig.  2.  Frequency  Response  of  Circularly  Symmetric  2-D  Filter. 
N = 21,  o»  = .50rr,  = - 60ir , Wt  = 111,  6=  0.055 


closely  approximate  the  true  optimal  filters  designed  by 
any  of  the  previous  methods.  Using  McClellan  Transforma- 
tions, filters  of  large  orders  can  be  designed  in  a matter 
of  seconds  on  a general  purpose  computer.  Filters  of 
127x12?  sample  points  have  reportedly  been  designed  (Ref 
10i240),  although  the  memory  requirements  would  be  consid- 
erable . 

Approach  of  this  Investigation 

This  investigation  develops  a computer-aided  design 
program  which  extends  the  capability  of  the  McClellan 
Transformation  to  design  2-D  FIR  filters  with  a single 
transition  band  of  a fairly  arbitrary  shape.  Included  with- 
in the  admittable  class  of  shapes  are  low  pass  and  high  pass 
filters  that  are  symmetric  about  both  the  and  &>£  axes , 
and  whose  single  transition  band  is  bounded  by  smooth,  mono- 
tonic curves  defined  by  the  functional  relationship  a»2  = F(cw^ 
in  the  region  [o , rr] x [o , tt]  of  the  «m^-6>2  plane.  Additionally, 
it  is  required  that  the  attenuation  and  transmission  bands 
be  of  constant  magnitude. 

The  algorithm  to  be  implemented  consists  of  four  major 
steps i 

1.  Define  the  2-D  frequency  response  by  specifying 

admittable  contours  in  the  plane  that  will  approx- 

imate the  desired  transition  band  of  the  2-D  filter. 

2.  Using  an  error  criteria,  determine  the  values  t^, 
t2,  t-j,  t^  that  will  produce  the  closest  mapping  of 
1-D  transition  band  frequencies  to  2-D  transition  band 
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contours . 

3*  Design  an  optimal  1-D  FIR  prototype,  with  an  appro- 
priate frequency  response,  that  will  be  transformed  into 
the  desired  2-D  filter. 

4.  Calculate  the  impulse  response  coefficients  and 
the  frequency  response  of  the  designed  2-D  filter. 

The  complete  details  of  the  method,  until  now  unavailable, 
are  logically  developed. 

The  Sequence  of  Problem  Development 

Chapter  II  reviews  the  McClellan-Parks  algorithm  for 
designing  optimal  LSI  FIR  1-D  digital  filters.  This  chapter 
is  included  for  two  reasons:  first,  the  Parks-McClellan 
algorithm  is  used  to  design  the  1-D  prototypes  required  for 
the  2-D  transformation;  second,  the  principles  in  the  algo- 
rithm serve  as  an  essential  introduction  to  some  of  the 
concepts  of  the  later  chapters.  In  Chapter  III  the  mathe- 
matics of  2-D  sequences  is  reviewed,  and  the  required  tools 
for  2-D  design  are  developed.  Although  methods  for  design- 
ing non-optimal  2-D  filters  exist,  the  optimal  2-D  filter  re- 
mains a stubborn  problem.  The  reasons  that  make  true 
optimal  filter  design  in  two  dimensions  so  difficult  are 
examined.  This  motivates  the  ultimate  selection  of  a sub- 
optimal  method  for  the  2-D  filter  design  algorithm. 

In  Chapter  IV  the  rationale  behind  the  McClellan  Trans- 
formation is  discussed,  and  the  nature  of  the  complex  mapping 
required  in  the  contour  design  is  examined.  It  is  shown  that 
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the  method  cannot  design  transition  band  contours  of  all 
shapes,  but  can  design  arbitrary  contours  within  a very 
useful  class.  To  produce  a non-trivial  mapping,  the  re- 
quirements for  constraints  on  the  elements  t-^,  t2»  t^,  t^ 
of  Eq  (1)  is  shown  necessary. 

A method  for  calculating  the  2-D  frequency  response 
and  the  2-D  impulse  response  coefficients  is  outlined  in 
Chapter  V.  Additionally,  the  optimality  of  filters  designed 
using  the  McClellan  Transformation  is  considered.  Chapter 
VI  illustrates  a few  typical  design  problems.  Chapter  VII 
draws  conclusions  on  the  investigation,  and  makes 
recommendations  for  future  related  research. 

Several  appendices  are  included.  Of  particular 
interest  are  Appendix  C,  a discussion  of  the  computer  pro- 
gram developed  in  this  investigation  and  Appendix  D,  a user's 
guide  to  this  program. 
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II.  1-D  Optimal  FIR  Filters 


This  chapter  begins  with  a brief  review  of  finite- 
extent  sequences  and  the  nature  of  a linear  phase  1-D  FIR 
filter.  Finite  duration  sequences  possess  certain  very 
useful  properties  for  filter  design.  Such  problems  as 
stability  and  implementation,  which  might  affect  HR  filters, 
do  not  arise  with  FIR  filters  (Ref  7 « 75-76)*  Additionally, 
FIR  filters  can  be  designed  to  have  exact  linear  phase.  This 
is  often  important  for  applications  where  frequency  disper- 
sion due  to  non-linear  phase  is  harmful  (for  example,  speech 
processing  and  data  transmission)  (Ref  7*76) . For  many 
years  no  efficient  method  for  the  design  of  1-D  optimal  FIR 
filters  was  available.  The  Parks-McClellan  algorithm,  to  be 
described,  effectively  solved  this  problem. 

To  best  understand  the  Chebychev  (optimal)  approxima- 
tion problem,  certain  necessary  mathematical  theorems  are 
reviewed.  It  is  then  shown  that  the  Parks-McClellan  algo- 
rithm is  a natural  extension  of  earlier  efforts  to  apply 
this  Chebychev  approximation  theory.  However,  the  resource- 
ful adaptation  of  the  efficient  2nd  Remez  Exchange  algorithm 
gives  the  Parks-McClellan  method  far  more  power  than  any 
previous  method. 

Finite-Extent  Sequences  and  Linear  Phase 

If  (h(n)>  is  a causal  finite  duration  sequence  defined 
over  the  interval  0<n<2N-l,  then  the  z- transform  of 
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{h(n)}  is 


2N-1 

HU)  = l h(n) z-n  (2) 

n=0 

and  the  Fourier  transform  of  (h(n)}  , with  period  2^,  is 

2N-1 

H(w)  = £ h(n)exp(-jwn)  (3) 

n=0 

The  condition  for  linear  phase, 

H(*>)  = G(co)exp(  j (A  + Bto) ) , G(w)  Real  (4) 

will  be  achieved  if  (Ref  7:77-78) 

h(n)  = h(2N  - 1 - n),  0<n<2N-l  (5) 

Parks  and  McClellan  showed  that  one  of  four  possible  cases 
for  linear  phase  is  a positive-symmetric,  odd-length  impulse 
response  with,  in  Eq  (4),  A = 0,  B = -(2N-l)/2.  (Ref  1:10). 
Because  the  term  exp  j (A  + Bw)  contributes  only  to  phase, 
{h(n)}  can  be  assumed  centered  at  the  origin  with  length 
N + 1 in  each  direction.  Fig.  3 illustrates  a typical  case. 

If  his  assumption  is  made,  the  frequency  response  can 
be  rewritten  as  a sum  of  cosine  functions  (Ref  7 » 81-82): 

N 

G(o»)  = h(0)  + 2 T h(n)cos<un  (6) 

n=l 

N 

= a(n)cosa,n  (7) 

n=0 
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where 


a(0)  = h(0) 

(8) 

a(n)  = 2h(n) , n ^ 0 

(9) 

Chebvchev  Approximation  Problem 

By  defining  a desired  frequency  response  D(w) , and  a 

0 

weighting  function  W(«),  then  the  weighted  approximation 
error  E («)  is 

E(«)  = W(w)[d(w)  - G(oi)  ] (10) 

The  Chebychev  approximation  problem  amounts  to  finding  the 
coefficients  |a(n)}  that  minimize  the  maximum  absolute 
error  E(w)  over  the  frequency  band  of  interest. 

To  properly  develop  the  Chebychev  approximation 
problem,  a certain  amount  01  mathematical  background  is  re- 
quired. It  is  useful  to  extend  the  well-known  concept  of  a 
polynomial  approximating  function  by  considering  some  linear 

combinations  of  prescribed  functions  g^,  g2, gn  continuous 

on  a fixed  space  X.  Their  linear  combinations  Zc^g^  are 
termed  generalized  polynomials  (Ref  11:?2).  The  Existence 
Theorem  guarantees  that  to  each  f e X there  exists  at  least 
one  generalized  polynomial  which  best  approximates  f (Ref 
U«20). 

For  certain  types  of  generalized  polynomials,  the 
characterization  of  the  best  approximation  can  be  stated 
very  conveniently.  The  functions  g1,g2 g_  are  said  to 
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satisfy  the  Haar  Condition  if  each  g1  is  continuous  and 
if  every  set  of  n vectors  of  the  form 

x = [g1(x),  g2(x) g^x)] 

is  independent  (Ref  11:45).  For  systems  satisfying  the 
Haar  Condition,  the  following  theorem  holds  (Ref  11:75) 


Alternation  Theorem.  Let  g, be  a 

system  of  elements  of  C [a.bf  satisfying  the 
Haar  condition,  and  let  X be  any  closed  sub- 
set of  [a,b].  In  order  that  a certain  gener- 
alized polynomial  P = 2c. g.  shall  be  a best 

approximation  on  X to  a given  f e C [xj  it  is 
necessary  and  sufficient  that  the  error  function 
r = f - P exhibit  on  X at  least  n + 1 alter- 
nations. Thus  r(x^)  = -r(x^_-^). 

This  theorem  is  very  important  in  numerical  determina- 
tion of  the  best  approximations.  To  illustrate  the  theorem, 
it  might  be  desired  to  approximate  the  function  sin('fx/2) 
by  a first  degree  polynomial  P(x)  = cQ  + c1x  over  [o,l] 

(Fig.  4).  By  the  theorem,  the  error  function  must  alter- 
nate at  least  three  times.  The  points  of  alternation  are 
0 and  1,  which  can  be  determined  by  inspection,  and  £, 
which  is  not  known.  This  situation  is  typical:  if  the 
location  of  the  extrema  of  the  error  function  were  known, 
the  approximating  polynomial  could  be  determined  by  solving 
some  linear  equations.  To  determine  the  location  of  the 
extrema,  however,  usually  involves  solving  some  non-linear 
equations  (Ref  11:75-76). 

If  Eq  (10)  is  applied  to  the  design  of  an  optimal 
bandpass  filter,  the  number  of  extrema  of  E («)  will  be  the 


(11) 
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sum  of  extrema  of  G(w)  plus  the  mandatory  extrema  associ- 
ated with  a constraint  at  each  band  edge  (Fig.  5) • In 
general,  the  extrema  at  the  band  edges  (except  at w = 0 and 
w = tt)  will  not  be  extrema  of  G(«)  (Ref  7:129).  "Extra- 
ripple" may  therefore  be  present,  since  more  than  n + 1 
alternations  might  exist. 

Development  of  an  Efficient  1-D  Design  Method 

By  specifying  N and  61  (,62  = k6j)  and  allowing  Uq  and 
a»p  to  be  free  variables,  Hermann  designed  extraripple 
optimal  FIR  filters.  He  was  restricted  to  designs  of  fairly 
low-order  filters  because  of  an  inefficient  method  for 
solving  the  non-linear  equations  which  locate  the  extremal 
frequencies  (Ref  2:258). 

Using  a method  of  polynomial  interpolation  "reminis- 
cent of  the  Remez  Exchange  algorithm",  Hofstetter  et  al 
provided  an  efficient  method  for  determining  these  extremal 
frequencies  (Ref  12).  This  avoided  the  restriction  on  the 
filter  order,  but  still  left  «s  and  cd^  as  free  variables. 

Parks  and  McClellan  completed  the  solution  by  develop- 
ing a method  where  N,  «s,  and  are  fixed,  and  6^  ( 62  = *6^) 
is  allowed  to  vary  (Ref  3»  Ref  4).  By  properly  interpreting 
the  Alternation  Theorem,  they  developed  an  algorithm  capable 
of  designing  optimal  FIR  filters  with  the  minimum  number  of 

i 

error  alternations.  Following  on  Hofstetter' s work,  they 
incorporated  the  very  efficient  2nd  Remez  Multiple  Exchange 
algorithm  (Ref  13)  to  provide  their  method  with  great  speed. 
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For  example,  a 100  sample  point  filter  can  be  designed  in 
about  20  seconds  (Ref  4:97).  If  6^  must  meet  some  tolerance, 
the  Parks-McClellan  method  can  be  repeated  with  adjustments 
in  one  of  the  parameters . 

The  source  listing  of  Parks  and  McClellan's  FORTRAN 
program  for  designing  optimal  FIR  filters,  as  well  as  an 
excellent  discussion  of  the  2nd  Remez  algorithm,  is  avail- 
able in  Rabiner  and  Gold's  text  (Ref  7:187-204).  With 
certain  minor  modification,  this  program  was  used  to  pro- 
vide the  1-D  prototypes  for  use  in  the  design  of  the 
transformed  2-D  filters. 

Summary 

The  Parks-McClellan  algorithm  provides  a very  fast 
method  for  designing  1-D  optimal  FIR  digital  filters.  Be- 
cause this  was  one  of  the  last  remaining  major  1-D  filter- 
ing problems,  it  is  not  surprising  that  research  in  the  past 
several  years  has  centered  on  higher-dimensional  filtering, 
especially  the  two-dimensional  case.  The  insights  gained  in 
considering  the  1-D  optimal  design  problem  prove  useful  in 
the  discussion  of  2-D  optimal  filter  design. 


r 
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III.  The  Design  of  Optimal  2-D  Filters 

Techniques  for  designing  1-D  digital  filters  are  now 
fairly  well  established  in  the  literature.  Very  few  of 
these  techniques,  however,  have  been  extended  to  two  or  more 
dimensions.  Since  there  is  no  lack  of  interest  in  2-D 
signal  processing,  it  is  worth  investigating  those  design 
techniques  that  can  be  extended  to  two-dimensions. 

Of  special  interest  in  this  investigation  is  the  design 
of  2-D  optimal  filters.  In  this  chapter,  the  design  methods 
currently  available  are  reviewed,  each  being  evaluated  in 
terms  of  its  optimality,  computational  efficiency,  and  ease 
of  understanding.  It  is  shown  that  no  optimal  method  meets 
all  three  criteria.  This  justifies  the  selection  of  a "sub- 
optimal”  design  method. 

Basic  Concepts  of  2-D  Discrete  Signals 

Because  1-D  LSI  systems  represent  a special  case  of 
multi-dimensional  systems,  it  is  not  surprising  that  many 
of  the  concepts  encountered  in  2-D  filtering  seem  familiar. 
However,  certain  useful  properties  are  unique  to  the  1-D 
case. 

Discrete  2-D  signals  are  functions  of  two  integer 
variables  (m,n).  A discrete  function  x(m,n)  * is  not  de- 
fined unless  both  m and  n assume  integer  values.  In 

* 

The  "hat"  symbol  is  used  to  maintain  consistency  with 
the  FORTRAN  program  in  Appendix  C,  where  two-dimensional 
arrays  are  identified  with  a "HAT"  suffix. 
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practice,  2-D  discrete  signals  will  often  be  sampled  values 
from  a 2-D  continuous  signal,  and  appropriate  2-D  sampling 
theorems  exist  (Ref  14) . 

As  in  1-D,  two  dimensional  LSI  systems  can  be  completely 

A 

specified  by  their  impulse  response  h(m,n).  The  output 
y(m,n)  of  a 2-D  LSI  system  is  the  convolution  of  the  input 
sequence  x(m,n)  with  h(m,n): 

y(m,n)  = £ £ x(k,f)h(m-k,  n-1 ) (12) 

k=-«o  |=-«o 

The  complex  exponentials  z™  and  z2  are  the  eigenfunctions 
of  LSI  systems.  Thus, 

y(m,n)  = z-j"1  z2n  H(z1,z2)  (13) 

A 

where  H(z^,z2)  is  called  the  system  function. 

a . . . 

The  z-transform  of  x(m,n)  is 

X(zltz2)  = £ £ x(m,n)z1_m  z2"n  (14) 

m=-«  n=-<» 

and  the  system  function  can  be  shown  to  be  the  z-transform 
of  the  impulse  response  (Ref  10«231).  If  the  z-transform 
is  evaluated  on  the  surface  z1  = exp(j«^),  z2  = exp(j<*»2), 
the  2-D  Fourier  Transform  can  be  evaluated 
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(15) 


x («,,«-•)  = £ 

• m=-«o 


oo 

Y x(m,n)exp(- j(o>,m  + to  n)) 

n=-« o x c 


which  is  periodic  over  [-Tr,Tr]  in  both  dimensions. 

As  in  1-D, 

«*>2)  = X(wltw2)  H (<«>1 , Wg ) (16) 

The  2-D  discrete  Fourier  Transform  (DFT)  is  defined  to  be 

a M-l  N-l  A 

X(k,  1)  = Y £ x(m,n)exp(- j(2rrmk/M  + 2rmi/N))  (17) 

m=0  n=0 

As  in  1-D,  the  2-D  DFT  has  proven  to  be  an  extremely  useful 
tool  in  implementing  digital  filters,  especially  when  an 
efficient  Fast  Fourier  Transform  (FFT)  is  used.  Although 
direct  convolution  can  be  used  to  implement  filters  of  low 
order,  it  becomes  worthwhile  to  use  an  FFT  for  orders  of 
roughly  10x10  sample  points  or  more  (Ref  15:405). 

The  2-D  Design  Objective 

The  normal  objective  of  2-D  filtering  is  to  alter  the 
frequency  spectrum  of  some  input  function  by  operating  on  it 

A 

with  a filter  function  H(«^,w2).  The  digital  filter  design 
problem  amounts  to  the  determination  of  the  impulse  response 
coefficients  h(m,n)  which  will  produce  a filter  capable  of 
performing  the  desired  operation.  Of  course,  after  the 
filter  is  designed,  it  still  remains  to  be  implemented. 
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2-D  HR  Filters . For  a 2-D  IIR  filter,  the  system 
function  will  take  the  form  of  a ratio  of  finite-degree 
polynomials  (Ref  10:233) 


A 

H(Z;l 


A ( z^ • z2 ) 
B ( Zl , ) 


= Z Z h(m,n)z1"m 

m=-«>  n=-« 


(18) 


where 

a M1  N1  ^ 

A(zltz2)  = X X a(k,i)  z,'k  z ~f  (19) 

^ k=0  1=0  x * 

M2  N2 

B(z,,z?)  = X Z Mk,i)  z “k  z~l  (20) 

k=0  1=0  ^ 


then  a recursive  partial  difference  equation  exists: 

M,  N, 

A . . . 1 A . . A . 

y(m-i,n- j)  = ryr — rv  Z Z a^k»  Ox(m-k.n-l) 

bU.j)  k“0  1^0 


m2  n2 

h(l  i)  Z Z b(k,l)y(m-k,n-l)  (21) 

ovi.j;  k=Q  1=0 
k/  1 =0 


The  order  of  recursion  in  Eq  (21)  is  important,  for  of  the 
many  possible  schemes  of  recursing  on  the  variables,  only  a 
few  will  be  causal.  As  a simple  illustration,  consider  the 
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following  partial  difference  equations 

y(m,n)  = by(m-l.n)  + x(m,n)  (22) 

where  several  possible  recursions  will  be  considered.  In 
Fig.  6 a few  discrete  points  near  the  origin  are  labeled 
and  the  boundary  condition 

A.  X 

y(m,n)  = 0,  m<0  ornsO  (23) 

is  indicated.  Causal  relations  would  be  described  by  any 
of  the  following  possible  recursive  paths: 

ABC ...  DEF ...  GHI ... 

ADG . • . BEH • . . CFI . . . 

A DB  GEC  ... 

However,  other  paths  can  be  selected  that  prevent  a causal 
relation: 

ABC. . • ■ • .FED  GHI. • . 

AEI ...  BF ...  DH ... 

AB  ED  GHI  FC  . . . 

A 2-D  HR  filter  can  be  implemented  from  a causal  partial 
difference  equation,  provided  the  required  boundary 
conditions  are  known  (Ref  16) . 

While  IIR  filters  offer  a much  richer  set  of  system 
functions  than  do  FIR  filters,  there  are  disadvantages 
which  make  recursive  filters  unattractive.  For  any  filter, 
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a prerequisite  for  stability  is  that  the  impulse  response 
be  bounded  (Ref  16) 


•o  oo  A 

£ £ h(m,n) 

Y*»— - _oq 


m=-oo  n=-«> 


FIR  filters  are  always  absolutely  summable,  due  to  their 
finite  extent,  so  they  are  always  stable.  HR  filters  must 
be  tested  for  stability. 


G oH  o1  o 


D oE  0F  o 


A /N  B 


/ S / S / / ' "i 

Boundary  Conditions  = 0 


Fig.  6.  An  Example  of  Possible  2-D  Recursions 
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In  1-D,  an  IIR  filter  of  the  form 


HU)  = |f|}  (25) 

can  be  tested  for  stability  by  factoring  the  functions  A(z) 
and  B(z)  into  cascades  of  first  and  second  order  terms.  The 
Fundamental  Theorem  of  Algebra  (Ref  17)  guarantees  that  this 
can  be  done.  Unfortunately  there  is  no  analogous  theorem 
for  polynomials  of  two  or  more  variables.  Stability  must 
be  tested  in  a more  indirect  way. 

The  most  practical  IIR  stability  test  is  due  to  Huang 
(Ref  18).  The  test  can  be  summarized  as  follows: 

' A 

Huang’s  Test.  A recursive  filter  H(z,,z?)  is  stable 
if  and  only  if  the  image  of  the  unit  circle  |Z2|  = 1 
when  mapped  into  the  z plane  by  6(zi,z2)  = 0 does  not 
intersect  the  region  |z-jj  > 1 and,  in  addition,  no 
point  in  the  region  lz,|>l  maps  into  the  point 
z2  = 0. 

Theoretically,  this  test  would  require  the  mapping  of  the 
entire  Z±-Z2  hyperspace.  In  practice,  only  a discrete 
number  of  points  can  be  mapped.  Improvements  on  Huang's 
Test  have  been  published,  but  the  test  remains  a major 
computational  problem  (Ref  10:235)* 

Even  if  a stability  test  is  available,  more  important 
problems  remain:  first,  how  the  designer  determines  the 
filter  coefficients  (of  a stable  filter)  to  meet  a pre- 
scribed frequency  response  specification;  second,  how  an  un- 
( stable  filter  can  be  stabilized  without  adversely  affecting 

magnitude  response. 
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Complex  cepstrum*  techriques  suffer  because  they  yield 
infinite  degree  denominator  polynomials  which  must  be 
approximated  by  finite  degree  polynomials,  sometimes 
altering  the  frequency  response  significantly  (Ref  10s235). 

A stabilization  conjecture  due  to  Shanks  (Ref  19) , that  the 
double  planar  least-squares  inverse  (PLSI)  of  an  unstable 
filter  would  yield  a stable  filter  with  the  desired 
characteristics,  has  been  used  to  design  many  HR  filters. 

For  years,  no  counter-example  appeared.  However,  very 
recently,  Kamp  and  Genin  provided  a counter-example  (Ref  20) 
and  then  disproved  the  conjecture  (Ref  21). 

2-D  FIR  Filters . Because  of  these  problems  with  2-D 
IIR  filters,  most  2-D  designs  have  used  FIR  filters.  If 
the  desired  2-D  filter  either  requires  or  can  tolerate 
circular  symmetry,  then  several  design  techniques  for  FIR 
filters  are  available.  By  circular  symmetry,  it  is  meant 
that  any  slice  of  the  filter  that  passes  thru  the  origin 
will  be  identical  (Fig.  7)* 

The  most  straight-forward  FIR  design  technique  is  the 
direct  extension  of  the  windowing  method  to  two-dimensions. 
Using  windows,  the  infinite-extent  2-D  Fourier  coefficient 
h(m,n)  are  multiplied  by  a finite-extent  2-D  window  function 
w(m  ,n)  to  produce  the  coefficients 

k(m,n)  = h(m,n)  *w(m,n) , -M<m<M,  -N  < n < N (26) 


Defined  as  the  inverse  Fourier  transform  of  the  complex 
logarithm  of  the  Fourier  transform  of  a sequence. 
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The  approximated  frequency  response  becomes  (Ref  16) 


a M N a 

H(o)  cj  ) = Y V k(m,n)exp(- j(w  m + w9n))  (27) 

1 m=-M  n=-N  x 

Huang  has  shown  that  "good"  2-D  windows  can  be  obtained 
from  "good"  1-D  windows  by  the  relation 

w(m,n)  = w (yj  m2  + n2  ) (28) 

where  w(»)  is  a continuous  1-D  window  (Ref  22).  Using 
windows,  circularly  symmetric  filters  with  a great  variety 
of  shapes  can  be  designed. 


Fig.  7.  Slices  from  a Circularly  Symmetric  Filter 
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There  are  drawbacks  to  windowing,  however.  First,  in 
the  design  of  a 1-D  prototype  window,  a closed  form 
expression  for  the  Fourier  series  coefficients  is  required: 

, 2tt 

h(n)  = ■£ -j  H(w)  exp  (jeon)  do  (29) 

If  this  integral  proves  cumbersome  to  evaluate,  the  design 
process  may  stop  (Ref  7*101).  Second,  window  designed 
filters  cannot  meet  prescribed  frequency  cutoff  specifica- 
tions because  the  convolution  process  "smears"  the  frequency 

response  near  the  ideal  response's  discontinuities.  Although 

{ 

there  may  exist  a "best"  window,  it  would  not  produce  an 
optimum  filter  (Ref  ?:90). 

Rabiner  and  others  have  used  linear  programming  to 
design  2-D  FIR  filters.  This  approach  can  be  considered 
a natural  extension  of  the  1-D  frequency  sampling  technique 
(Ref  7*105-123) • except  that  all  the  frequency  samples  are 

A 

considered  as  variables.  If  Dtw-^.wg)  is  the  desired  ideal 
response,  then  the  constraints 

A A 

H(g»i,<*>2)  < EK^i'Ci^)  + 6 (30) 

A A 

H (a>i , g>2 ) > D(o)1,o)2)  - 6 (31) 

can  be  applied  to  a grid  of  points  in  the  P^3116* 

A 

Linear  programming  adjusts  h(m,n)  while  minimizing  6 . 

Optimal  2-D  filters  have  been  designed  using  linear  pro- 
gramming, but  the  method  proved  computationally  costly. 
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Rabiner  reports  that  the  most  ambitious  design  he  attempted 
(a  9x9  sample  point  filter  with  circular  symmetry)  required 
more  than  one  hour  of  execution  time  on  a high  speed  com- 
puter (Ref  7:461-471) . Since  thousands  of  constraint 
equations  can  be  involved  in  the  design  of  a 2-D  optimal 
filter,  linear  programming  (which  is  a single -exchange 
algorithm)  would  be  expected  to  be  an  inefficient  design 
method. 

Efforts  have  been  made  to  extend  the  design  method  of 
Parks  and  McClellan  to  two-dimensions.  McClellan  considered 
the  problem,  but  stopped  after  realizing  that  the  2nd  Remez 
algorithm  could  not  be  extended  to  the  2-D  case.  Recently, 
in  independent  efforts,  Kamp  and  Thiran  (Ref  8)  and  Kersey 
and  Mercereau  (Ref  9)  developed  2-D  design  techniques 
similar  to  the  Parks -McClellan  approach  but  based  on  the 
1st  Remez  Multiple  Exchange  algorithm  (Ref  13).  Since  a 
multiple  exchange  algorithm  was  used,  an  impressive  speedup 
in  comparison  with  linear  programming  was  hoped  for.  How- 
ever, the  speedup  was  not  dramatic,  because  the  bulk  of  the 
execution  time  was  spent  in  the  evaluation  of  the  error 
function,  rather  than  in  optimization  of  parameters  (Ref 
10«240) . 

Analysis  of  these  methods  reveal  that  the  Chebychev 
approximation  of  2-D  FIR  filters  is  fundamentally  more 
difficult  and  slower  than  the  1-D  case.  There  are  three 
reasons.  First,  there  is  the  impossibility  that  any  set  of 
functions  defined  on  a 2-D  domain  can  satisfy  the  Haar 


condition,  thus  weakening  the  Alternation  Theorem  consider- 
ably (Ref  9*8).  Second,  it  would  not  be  possible  to  order 
the  extremal  frequencies,  as  in  1-D,  where  increasing  order- 
ing guarantees  that  the  error  sign  alternates  from  point  to 
point.  Third,  the  size  of  the  problem  would  approach  the 
unmanageable  as  the  number  of  optimization  parameters  is 
increased  (Ref  1:30).  In  the  absence  of  the  Haar  condition, 
the  1st  Remez  algorithm  must  be  used,  and  convergence  may 
not  occur  unless  suitable  perterbations  are  introduced 
(Ref  8:262).  Lack  of  a simple  alternation  in  the  error 
complicates  the  rules  for  exchanges,  at  times  preventing 
multiple  exchange,  the  very  procedure  that  makes  the  Remez 
algorithms  efficient. 

As  with  all  the  design  techniques  presented  so  far, 

Kamp  and  Thiran  limited  their  designs  to  circular  symmetry. 
Hersey  and  Mercereau  improved  on  this  by  requiring  only 
octagonal  symmetry*.  Although  they  were  able  to  design 
filters  with  a more  arbitrary  frequency  response,  computa- 
tional considerations  limited  their  filters  to  15x15  sample 
points  or  fewer  (Ref  9) . 

Because  these  methods  design  truly  optimal  filters, 
they  have  definite  value.  However,  their  costly  execution 
time  confines  their  use  to  relatively  small  impulse  responses. 
This,  coupled  with  their  mathematical  sophistication, 
probably  prevents  them  from  being  widely  accepted. 


Symmetry  about  both  axes  and  the  line  «»2  = 


Summary 


Until  optimal  methods  become  computationally  competitive 
with  non-optimal  methods,  most  2-D  filters  will  probably 
be  designed  non-optimally.  Fortunately,  a new  method  pro- 
posed by  McClellan  is  capable  of  improving  significantly  on 
previous  non-optimal  methods.  It  designs  filters  that  are 
equiripple,  although  not  quite  optimal.  As  such,  they  have 
been  called  "sub-optimal”  (Ref  9*1) • This  new  method  is 
discussed  in  the  next  chapter. 


( 
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IV . The  McClellan  Transformation 


A new  method  for  designing  2-D  digital  filters,  due 
to  McClellan  (Ref  1),  has  been  cited  as  the  design  method 
of  choice  for  most  applications  today  (Ref  9:1-2). 

The  essence  of  the  method  is  that  it  transforms  an  optimal 
1-D  FIR  filter  into  a sub-optimal  2-D  FIR  filter.  These 
filters  are  termed  "sub-optimal"  because  they  retain  the 
equiripple  characteristic  of  true  optimal  filters,  yet  they 
are  inferior  to  optimal  filters  when  judged  very  strictly. 

As  mentioned  in  Chapter  I,  McClellan  has  designed  circularly 
symmetric  bandpass  filters.  However,  he  proposed  that  his 
method  could  be  used  to  design  filters  with  a much  less 
restrictive  symmetry  (Ref  1:42-43).  The  important  features 
of  the  method  are  that  it  is  not  limited  to  the  design  of 
small  filters  (127x12?  point  filters  have  been  designed) ; 
that  it  is  fast,  efficient,  and  simple  to  understands  and, 
surprisingly,  in  many  cases  the  resulting  filters  are  almost 
indistinguishable  from  true  optimal  filters  (Ref  24:93-97). 
This  last  point  will  be  considered  again  in  detail  in 
Chapter  V. 

Spectral  Transformation  Mappings 

McClellan's  method  can  be  considered  a spectral 

transformation  mapping.  Mitra  provides  this  definition 

of  a spectral  transformation  (Ref  25*905)* 

By  a spectral  transformation,  we  mean  a complex  map 
that  carries  a stable  rational  transfer  function 
into  another  stable  rational  transfer  function 


exhibiting  a different  frequency  response,  at  the 
same  time  maintaining  some  desirable  characteris- 
tics... In  general,  a spectral  transformation  must 
have  the  following  characteristics:  1)  it  must 
produce  a LSI  stable  transfer  function  from  an  LSI 
stable  transfer  function;  2)  it  must  transform  a 
real  transfer  function  into  a real  transfer  function; 
and  3)  it  must  preserve  some  basic  characteris- 
tics of  the  magnitude  response  (e.g.  the  ripple 
magnitude  in  the  transmission  and  attenuation 
domain) . 


The  spectral  transformation  that  McClellan  adopted  was 
first  stated  in  Chapter  I and  is  repeated  here: 

COSU  = t^  + t2COSw^  + t^COSu^  + t^COStJ^*  COSW2  (1) 

where  u,  and  w2  are  restricted  to  the  interval  [-Tr,Tr], 
as  discussed  in  Chapter  II.  Eq  (1)  is  referred  to  as  the 
McClellan  Transformation. 

The  Change  of  Variables 

In  Chapter  II  it  was  shown  that  a 1-D  linear  phase 
filter  of  length  2N  x 1 has  the  form 

H(«)  = G(w)exp( j (A  + Bw))  (32) 

where  the  term  exp  j (A  + Bcu  ) contributes  only  to  phase 
and  can  be  disregarded.  (Otherwise  it  would  be  carried  along 
unaltered  through  each  transformation  to  be  described.)  It 
followed  that 

N 

G(w)  = £ a(n)coswn  (33) 

n=0 

c 
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The  change  of  variables  in  Eq  (1)  depends  on  G(«)  being 
expressable  by  a suitable  choice  of  coefficients  {b(n)}  as 
the  following  linear  combination  of  cosine  functions: 

N N 

G(eo)  = T,  a(n)coswn  = £ b(n)*(cosw)n  (34) 

n=0  n=0 

4 

To  calculate  fb(n)}  , let  x = cosw.  Then 

Tn(x)  = cos£n  cos-^(x)J  = cosa>n  (35) 

where  Tn(x)  is  the  nth  order  Chebychev  polynomial  of  the 
first  kind.  Thus, 

N N N 

G(o>)  = V a(n)cosoin  = T a(n)T  (x)  = V b(n)xn  (36) 

n=0  n=0  n n=0 

and  by  expanding  each  equation  in  Eq  (36) 

H(«)  = a(0)TQ(x)  + a(l)T1(x)  + ...  + a(n)Tn(x) 

= b(0)x°  + bdJx1  + ...  + b(n)xn  (37) 

and  then  using  the  Chebychev  recursion  formula 

Tn+1(x)  = 2xTn(x)  - Tn-1(x)  (38) 

one  can  solve  for  x°,  x^-,  x2,...,xn  (the  method  to  do  this 


is  described  in  Appendix  A) « 


( 


or 
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(39) 


x = u T (40) 

0 1 

By  substituting  for  x , x ,...,xn  and  grouping  common 
factors  of  Tq(x),  T-^(x) , . . . ,Tn(x)  in  the  right  half  of  Eq  (37), 
one  can  equate  coefficients  of  Tq(x),  T^(x) , . . . ,T  (x)  to 
yield 


a(0) 

. a(l) 

•a(2) 

a(3) 

• 

= 

• 

• 

a(n) 

0.5  0 

1 0 .75 

0.5  0 

0 0 .25 


or 


Thus 


A = UT  B 


B = (UV1  A 


b(0) 

b(l) 

b(2) 

b(3) 

• 

b(n) 


(41) 


(42) 


(43) 
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which  allows  b(0),  b(l),...,  b(n)  to  be  solved  whenever  U 
is  non-singular. 

By  substituting  Eq  (1)  into  the  very  convenient  form 
of  Eq  (34), 

a II  r 

G(&>jL,w2)  = 2,  b(n)  t1  + t2cos<«>1  + t^cosa^ 

+ tjjcoso^cosa^j11  (44) 

and  completing  the  expansion,  G(u>^,u>2)  can  be  expressed  in 
this  form: 

a N N A 

G(ci>^,a>2)  = £ £ a(m,n)  (cosu^m)  (cosw2n)  (45) 

m=0  n=0 


where  the  impulse  response  coefficients  are  directly  related 
to  a(m,n)  . The  methods  that  take  Eq  (44)  into  the  form  of 
Eq  (45)  are  involved  and  will  be  described  in  Chapter  V. 

It  is  noted  that  Eq  (45)  is  the  desired  form  for  a 2-D 
FIR  zero-phase  digital  filter  (Ref  15:406).  This  is  what 
makes  McClellan's  particular  spectral  transformation  so 
useful.  In  fact,  it  is  possible  to  show  that  there  is  an 
entire  set  of  suitable  spectral  transformations 


P Q A 

cos«rt  = V T t(p,q)  (cospw,  ) (cosqaO 
p=0  q=0  1 * 


of  which  Eq  (1)  is  the  special  case  P = 1,  Q = 1.  Eq  (46) 
has  been  called  the  "generalized  McClellan  Transformation" 
(Ref  24 » 20-24).  To  limit  the  scope  of  this  investigation, 
only  the  case  P = 1,  Q = 1 has  been  considered. 


The  Range  of  the  Mapping 

The  McClellan  Transformation  defines  a mapping  of  the 
interval  [-Tr,Tr]  of  the  1-D  frequency  axis  to  the  square 
region  [-tt,tt]x  [-Tr.Tf]  of  the  2-D  frequency  plane.  The  map- 
ping is  certainly  not  one  to  one , and  the  range  of  the 
mapping  may  not  include  the  entire  square  region  in  . 

If  quadralateral  symmetry* is  imposed,  the  problem  can  be 
equivalently  stated  as  the  mapping  [O , tt]  — [o , tt]  x[0,tt], 
since  the  information  in  the  first  quadrant  of 
repeated  in  the  other  three  quadrants. 

If  «2  is  a one-valued  function  of  then  Eq  (1)  can 


be  solved  for  u»2: 


a>2  = cos 


COS  Cl)  — t-j^  t2COSCJ11 

|_  + t^COSo»2  i 


For  each  fixed  value  of  c []o,tt]  there  exists  some 
corresponding  contour  C^  in  the  «^-w2  Along  this 

contour,  the  magnitude  of  the  2-D  frequency  response  is 
identical  to  the  1-D  frequency  response  at  By  allowing 

a>  to  vary  from  0 to  tt,  a complete  family  of  contours  is 


* 

Defined  to  be  symmetry  across  the  and  a>2  axes. 
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generated  describing  the  2-D  frequency  response  over  some 
portion  of  the  region  [0 , tt]  x [0 , tt  J . For  example,  in  the 
original  McClellan  transformation,  the  values  t^  = -.5, 
tg  = tj  = t^  = 0.5  produced  the  circularly  symmetric 
contours  shown  in  Fig.  1.  For  frequencies  below  0.7tt, 
the  contours  are  very  close  approximations  to  circles.  How- 
ever, as  cd  — tt,  the  contours  become  more  squarelike. 

This  is  still  a very  desirable  mapping,  because  it 
maintains  the  desired  contour  shape  over  a large  area  of  the 
2-D  region,  and  yet  there  is  a correspondence  between  every 
point  in  £o,ttJx  [0,1/]  and  the  1-D  interval  [o.ttJ.  Actually, 
forcing  the  mapping  to  have  this  range  will  almost  always 
force  some  distortion  in  the  contour  shape.  But,  if  the 
mapping  is  not  so  defined,  the  resulting  frequency  response 
of  the  2-D  filter  will  be  ill-behaved.  This  is  discussed 
in  detail  later  in  this  chapter. 

The  Design  Strategy 

The  design  strategy  and  its  limitations  are  now  fairly 
clear.  A 1-D  filter  such  as  shown  in  Fig.  5»  with  transi- 
tion band  frequencies  cd^  and  &>s,is  mapped  into  a 2-D  filter 

filter  with  transition  contours  C and  C . The  magnitude 

p s 

response  in  the  transmission  and  attenuation  domains  is 
carried  unaltered  from  one  dimension  to  two  dimensions.  If 
the  1-D  prototype  filter  was  an  equiripple  optimal  filter,  then 
the  2-D  transformed  filter  will  also  be  equiripple.  How- 
ever, it  will  only  be  optimal  along  the  slice  cd2  = 0. 


38 


( 


Everywhere  else  it  will  be  slightly  sub-optimal  (Ref  24s 
93-97).  Because  this  is  usually  a small  trade-off  in 
comparison  with  the  relative  ease  and  flexibility  of  design, 
the  McClellan  transformation  is  a very  attractive  method 
for  designing  2-D  FIR  filters. 

There  are  limitations  on  the  types  of  contours  that 
can  be  mapped  using  the  method.  This  can  be  seen  by  letting 
x = cos  <i>  , u = cos  and  v = cos  a>2,  then  Eq  (47)  becomes 


x - t1  - t2u 
t3  + t^u 


(48) 


If  x is  held  fixed  and  a partial  differentiation  is  per- 
formed (Ref  1:36-38) 


av  = ~t2t3  + tlt4  1 V (4 

9U  (t3  + t4u)2 

it  is  seen  that  v is  monotonic,  since  the  sign  of  the 
derivative  does  not  change  as  u varies  from  -1  to  1.  This 
requires  that  the  contours  in  a>^-«2  must  also  be  monotonic. 

The  Contour  Approximation  Problem 

In  order  to  design  2-D  filters  with  arbitrary  contour 
shapes  within  the  admissable  class,  some  method  must  be 
developed  that  maps  prescribed  1-D  frequencies  to  a speci- 
' fied  2-D  contour.  As  a first  cut  at  the  problem,  it  is 
assumed  that  only  one  of  the  1-D  transition  band  frequencies 
must  be  mapped  to  a particular  contour.  Hopefully,  the 
contour  shapes  in  the  neighborhood  of  the  designated 
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contour  will  not  be  greatly  changed.  Let  the  1-D  frequency 

be  , and  the  desired  contour  be  C , where  C is  defined 
o o o 

by  the  single  valued  relation 

"2  = Fo(<ul)  (5°) 

Then  for  some  value  of  t 

t = 


o>0  must  be  made  to  map  to  Cq.  In  general,  this  is  impossi- 
ble, so  a best  approximation  according  to  some  error  criteria 
is  selected.  Two  of  the  possible  criteria  are  to  minimize 
the  maximum  absolute  error  (minimax  criteria)  and  to  mini- 
mize .c  sum  of  the  squared  errors  (least-squares  criteria).  ' 
Because  of  the  relative  ease  of  implementation,  the  least- 
squares  criteria  was  selected.  Empirical  evidence  indicates 
that  the  two  methods  lead  to  equivalent  filters,  at  least 
for  the  case  P = 1,  Q = 1 in  Eq  (46)  (Ref  2«407-408). 

By  considering  the  constant  cosw0  to  represent  a set 
of  identical  sample  values,  the  classical  least-squares 
formulation  for  fitting  a set  of  data  points  can  be  used. 

Let 

S^)  = t1G1(oi1)  + t2G2(w1)  + t^G^o^)  + t^G^^)  (52) 


(5D 
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where 


= cosa>o  = constant  (53) 

G1(o,1)  = 1.0  (54) 

G2(wi)  = cos^  (55) 

G'3 (t«>i ) = cos(  F(a»1)  ) (56) 

^4 (W1 ) = c05^!)  cos(  F0(ft'l)  ^ (57) 


A requirement  that  each  of  the  approximating  functions  (in 
this  case,  they  are  generalized  polynomials)  be  linearly 
independent  is  satisfied.  The  least-squares  formulation 
provides  four  equations  to  solve  for  the  four  unknowns  t. 

Because  an  explanation  of  the  least-squares  formulation 
is  lengthy  and  might  be  unnecessary  for  many  readers , the 
details  of  this  method  as  applied  to  the  contour  approxi- 
mation problem  are  outlined  in  Appendix  B.  References  on  the 
method  are  included  in  the  bibliography  (Ref  27,  Ref  28). 

Constraints  on  the  Mapping 

This  much  accomplished,  the  contour  approximation 
problem  still  remains  only  half  solved.  For  example,  if  t 
is  free  to  take  on  all  values,  the  trivial  solution 
t-j_  = cosa>o,  t2  = t^  = t^  = 0 will  result  in  a zero  error. 

> This  would  also  result  in  the  frequency  &>o  mapping  to  the 
entire  2-D  region  [0 , tt] x [0 , tt]  as  well  as  CQ.  Evidently, 
some  reasonable  constraints  must  be  placed  on  t to  reach  a 
useful  solution. 
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A possible  set  of  constraints  that  seems  intuitively 
satisfactory  for  the  mapping  — [0 , tt]  x [0 , tt]  might  he 


(0,0) 

(58) 

(tr , Tf ) 

(59) 

G(co2,o)1) 

(60) 

tr  — 

A 

G (<«>i , o>2 ) = 

where  the  second  two  constraints  imply  a redundant  constraint 

tt  — (tt,  2)  (61) 

Applying  the  first  constraint  to  the  McClellan  Transforma- 
tion yields 


1 - t^  + t2  + t^  + t^ 


while  the  second  constraint  yields 


(62) 


— 1 — t^  — t^  + (t2  — t^)  cosa>-^ 

Because  Eq  (63)  must  he  true  for  every  e [O.ttJ,  it  is 
also  true  that 


(63) 


-1  = tx  - t3 

(64) 

0 = t2  - t4 

(65) 

Finally,  the  constraint  of  Eq  (60)  requires  that 


r 


t2  = t3 


(66) 


as  seen  hy  inspection  of  the  symmetry  of  Eq  (1) . Solving 
Eqs  (62) -(66)  for  t yields  = -.5,  t£  = t^  = ti+=  0.5, 
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which  is  recognized  as  the  original  McClellan  transformation 
for  circular  symmetry. 

A problem  with  this  set  of  constraints  is  that  it 
leaves  no  free  variables.  To  design  arbitrary  contours,  a 
less  restrictive  set  of  constraints  is  needed.  Mitra's 
paper  on  spectral  transformations  provides  two  theorems 
that  are  helpful  (Ref  25*910),  but  a superior  technique  has 
been  introduced  by  Mercereau  et  al  (Ref  15*408-409)  and  is 
recommended  for  most  contour  designs.  Their  suggestion 
is  to  initially  perform  the  contour  approximation  subject 
to  a single  constraint,  0 — (0,0).  This  produces  a map- 
ping to  an  intermediate  domain  (n^,^)-  A second  mapping 
will  produce  the  final  result*  a mapping  in  the  (cj-^.a^) 
plane  that  approximates  the  desired  contour  CQ  and  has  a 
range  that  includes  the  entire  region  [o , tT|  x [o , tt]  . 

As  an  example,  if  the  elliptical  contour 

o>2  = “>b  J 1 ■ (Wj/f^)2  (67) 

where 

(68) 

(69) 

(70) 


= -50  TT 

<“a  = .25* 

V 0 = *50x1 


( 
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is  desired,  the  least-square  formulation  subject  to 
constraint  (explained  in  Appendix  B)  yields 


-L  = -2.584 

(71) 

;2  = 2.586 

(72) 

3 = 0.166 

(73) 

;4  = 0.832 

(74) 

with  a mean  error  of  0.0003.  However,  the  range  of  the 
contour  mapping  is  inadequate,  as  shown  in  Fig.  8.  The 
significance  of  the  cross-hatched  area  can  be  understood  by 
considering  Eq  (1).  To  stay  within  the  realm  of  real 
numbers,  the  following  double  inequality  must  hold: 

-1.0  < costa  < 1.0  (75) 

or  equivalently 

-1.0  < < 1.0  (76) 

-l.o  < t-^  + t2cosn1  + t^cosn2  + t^cosn^cosr^  < 1.0  (77) 

The  cross-hatched  area  of  Fig.  8 corresponds  to  values  of 
( w^ere  Eq.  (77)  does  not  hold.  The  magnitude  of  the 

frequency  response  grows  rapidly,  since  it  corresponds  to 
complex  values  of  o>  (Refl5«408). 

Two  of  the  basic  operations  of  complex  variable  mapping 
theory  are  translation  and  expansion  (or  contraction). 

These  operations  do  not  alter  the  shape  of  a given  contour, 
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only  its  orientation  and  relative  size  (Ref  17s73)»  The 
following  mapping  should  be  considered* 

• At  A 

cosg>0  = F i>2)  = k1F(S21,n2)  - k2,  k-^  > 0 (78) 

The  usefulness  of  this  second  mapping  stems  from  considera- 
tion of  Eq  (77)  and  Fig.  8.  Because  the  filter  is  either 
lowpass  or  highpass , the  only  contours  of  importance  are  CQ 
and  its  immediate  neighbors,  since  they  control  the  mapping 
of  the  transition  band.  It  is  desired  that  the  range  of  the 
second  mapping  include  the  entire  square  region,  for  then 
Eq  (77)  would  be  satisfied.  Due  to  the  nature  of  the  second 
mapping,  the  contour  shapes  themselves  will  be  unaffected; 
however,  a new  1-D  frequency  will  be  associated  with  each 
contour  in  the  ultimate  mapping  domain. 

To  find  the  necessary  values  of  k^  and  k2  in  Eq  (78), 
one  applies  the  second  mapping  to  Eq  (76)  and  specifies 
the  constraint  of  Eq  (75): 

-1.0  < k-jF^,^)  - k2  < 1.0  (79) 

Solving  each  inequality  separately  for  its  worst  case  pro- 
duces two  equations  in  k^  and  k2: 


-1.0  - " k2 

(80) 

1.0  = k1Ffnax(n1.n2^  " k2 

(81) 

where  Fmin^itfi2)  and  Fmax(n1  ,n2)  correspond  to  the  least 


and  greatest  magnitudes  taken  by  F(  S2^,  n2)  over  a dense 
grid  of  points  in  £0 , tt]  x [0 , ttJ  . By  solving  Eq  (80)  and 
Eq  (81)  one  finds  that 


Fmax(“l,n2}  " Fmin(JVfi2^ 


Vmax^l’^V  ‘ 1 


(82) 

(83) 


The  new  1-D  frequency  associated  with  Cq  is  found  by 
solving  Eq  (78)  : 


too 

• .1 

= cos  (k-^cosa>0  - k2) 

(84) 

Because 

F (c^ , a>2 ) = 

• t • • 

ti  + t2  COSt^  + t^  COS&>2  + t^  COSto^ 

COS(U2  (85) 

a new  "scaled" 

t'  also  results: 

tl  = kltl  " k2 

(86) 

V = klt2 

(87) 

t3  = Vs 

(88) 

V " klt4 

(89) 

Using  the  t ' , a new  contour  mapping  is  shown  in  Fig.  9. 
It  is  noted  that  the  shapes  of  the  contours  remain  unchanged 
up  to  the  original  contour  C^.  In  o>1-o)2  the  contour 
now  corresponds  to 


Wl-flXIS  Cx  TO 


Fig.  9.  Elliptical  Contours  of  Eq  (67)  after  scaling. 

u0'  = 0. 2498-rr,  tT  = [-.0482,  .7559,  .0482,  .244l], 
i = .0002,  constraints  0^(0,0)  with  scaling 


(90) 


a»o  = cos"1  (kj^cosrr  - k2)  = cos"1  (-k^  -k2) 


The  remaining  contours  spread  out  over  fo.rr^x  [0,tt],  giving 
the  final  mapping  the  desired  range.  Thus,  the  final 
frequency  response  is  well-defined. 

As  mentioned,  several  theorems  attributed  to  San jit 
Mitra  can  also  be  helpful  in  assuring  a well-defined 
frequency  response.  Because  the  development  of  these 
theorems  is  extensive,  interested  readers  are  referred  to 
Mitra's  paper  (Ref  25)*  In  many  cases,  the  resulting  con- 
straints will  produce  in  the  same  final  mapping  as  in  the 
previous  section,  but  this  is  not  always  so.  To  leave 
the  design  program  as  flexible  as  possible,  either  method 
of  specifying  constraints  is  permitted.  The  technique  of 
Mercereau  et  al  is  available  as  a default  option;  however, 
an  enlightened  user  may  override  the  default  and  supply  any 
consistent  set  of  constraints. 

Summary 

The  McClellan  Transformation  is  a spectral  transforma- 
tion that  maps  1-D  frequencies  in  the  interval  to 

contours  in  the  (a»^,a>2) -plane  over  [0,TrJx[0,Tr]|.  Because  the 
filter  is  assumed  to  have  quadrilateral  symmetry,  the  area 
£ O.TfJxQO.Tr]  actually  defines  the  total  behavior  over 
Q-rr ,rr]  x [-tr,TfJ.  Because  the  mapping  is  generally  imperfect, 
an  error  criteria  is  established  to  determine  the  best 
mapping.  To  insure  that  the  mapping  is  non-trivial , con- 
straints on  t are  necessary.  To  insure  that  the  range  of 
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the  mapping  includes  the  entire  region  [o.ttJ  x Jo.ttJ  , the 
technique  of  Mercereau  et  al  can  he  used.  If  the  equiva- 
lence of  Eqs(44)  and  (45)  can  be  established,  the  filter 
design  is  complete. 


( 
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V.  Determination  of  the  2-D  Impulse 
Response  Coefficients 

As  shown  in  the  last  chapter,  t can  be  found  that 
provides  a suitable  contour  mapping.  Thus,  the  frequency- 
response  can  be  evaluated  as  a function  of 

a N 

= £ b(n)  (t1+t2costo1+t^cosw2+t^cosw^cosa)2)n  (44) 

where  w2  = F(cd^).  This  itself  is  not  enough,  since  it  is 
necessary  to  convert  Eq  (44)  into  a form  that  relates 
directly  to  the  impulse  response  coefficients  (h(m,n)}  . 

A 

They  will  be  needed  if  the  z-transform  K(z^,z2)  will  be 
used  when  the  filter  is  implemented. 

In  this  chapter,  a method  is  derived  that  transforms 
Eq  (44)  into  the  form  of  a zero-phase  2-D  FIR  filter.  From 
this  farm,  the  impulse  coefficients  can  be  determined  by 
inspection.  The  chapter  concludes  with  some  comments  on 
the  "optimality"  of  the  McClellan  Transformation. 

Algorithm  for  the  Expansion 

Let  x = coso>^  and  y = coso>2,  then 

A N „ N 

H(x,y)  = ^b(n)(t1+t2x+t3y+t4xy)  = £ b(n)[*]n  (91) 

The  evaluation  of  Eq  (91)  can  be  accomplished  by  consider- 
( ing  the  following  steps i 
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H0(x,y) 

= b(0) 

(92) 

H-^x.y) 

= Ho(x,y)  + 

b(l)[«] 

(93) 

H2(x,y) 

= H^(x,y)  + 

(94) 

H3(x,y) 

• 

A 

= H2(x,y)  + 

(V 

1 — 1 

• 

1 — 1 

1 — 1 

• 

1 1 

r>i 

(95) 

• 

• 

Hk(x,y) 

A 

= Hk_1(x,y) 

+ 

o' 

P? 

i — 1 

• 

i — i 

• 

i — i 

tv 

1 

H- < 

36) 

Because  the  quantity  ^t-^  + t2x  + t^y  + is  already 

known,  the  process  is  recursive.  As  the  order  of  N gets 
large,  the  memory  required  for  the  expansion  grows  very 
quickly.  However,  if,  after  each  step  of  the  expansion, 
all  of  the  coefficients  are  searched  for  like  terms , the 
memory  requirements  are  held  down  considerably.  Empir- 
ically, it  was  found  that  the  expression 

[tx  + t2x  + t3y  + t^xy]  [tx  + t2x  + t3y  + t^xy  ] k_1  (97) 

would  require  (k  + 1)  words  of  memory  for  coefficient 
storage,  provided  like  terms  were  combined.  If  the  search 
for  like  terms  was  delayed  until  the  (N  + l)th  step,  the 
storage  required  would  grow  as  4 . Table  I compares  the 
storage  required  with  the  two  possible  approaches. 

Obviously,  searching  for  like  terms  after  each  step  of  the 
expansion  is  essential  except  for  very  low  order  filter 
designs.  The  final  result  of  the  expansion  is 
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A N N A 

H(«,  ,<a_)  = £ Y b(m,n)  (cosca,  ) (cosu?) 

■L  c m=0  n=0 


(98) 


which  involves  every  possible  combination  of  (cos<a-^)m  (costa,,  )n 
up  to  m = N,  n = N. 

An  Analogous  1-D  Result 

Eq  (98)  itself  still  does  not  present  a direct  method 
of  obtaining  the  impulse  response  coefficients.  Before 
progressing,  it  is  instructive  to  consider  the  analogous 
1-D  equation 

a N _ N 

H(w)  = y b(n)  ( co sea;  = £ a(n)coscan  (99) 

n=0  n=0 


where  the  coefficients  a(n)  must  be  determined.  Using 
the  dual  of  the  argument  that  led  to  Eq  (43)  in  Chapter  IV, 
let  x = costa  then  Tn(x)  = coscan,  and 


Table  I 


STORAGE  REQUIRED  FOR 

TWO  METHODS 

OF  EXPANDING  EQ  (91) 

k 

(k  + l)2 

4k 

0 

1 

1 

1 

4 

4 

2 

9 

16 

3 

16 

64 

• 

• 

• 

• 

• 

• 

• 

• 

• 

10 

121 

1,048,576 

• 

9 

• 

• 

9 

• 

• 

9 

• 
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i 

(100) 


N „ N „ N 

H(«)  = Y b(n)  fcosw)  = Y b(n)xn  = Y a(n)T  (x) 
n=0  n=0  n=0  n 


Using  Eq  (38),  the  relationship  between  TQ(x),  T-^(x) 
0 1 n 

Tn(x)  and  x , x x can  be  shown  to  be 


or 


T0(x) 

l^fx) 

T2(x) 

= 

T3(x) 

• 

• 

• 

Tn(x) 

1 

0 

-1 

0 


0 

1 

0 

-3 


0 

0 

2 

0 


0 

0 

0 

4 


n 


(101) 


T = V X 


(102) 


( 


By  substituting  for  TQ(x),  T1(x) , . . . ,Tn(x)  and  grouping 

A 1 

common  factors  of  x , x , . . . ,x  in  the  right  half  of 
Eq  (100),  one  can  equate  the  coefficients  of  x°,  x^,...,xn 
to  yield 


b(0) 

b(l) 

b(2) 

= 

b(3) 

• 

• 

• 

b(n) 

0-10 
1 0 -3  ... 

0 2 0 

0 0 4 


a(0) 

a(l) 

a(2) 

a(3) 

• 

• • 

a(n) 


(103) 
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or 


f 


B = VT  A 


(104) 


Thus 


A = (VT)-1  B 


(105) 


which  allows  (a(n)}  to  be  solved  whenever  VT  is  non- 
singular. 


Extension  to  Two-Dimensions 

This  1-D  method  for  calculating  |a(n)}  can  be  exploited 
to  calculate  |a(m,n)}  . Eq  (98)  is  rewritten  as 


m n 


* N N A 

H(x,y)  = £ Y Mm,n)x  y' 

m=0  n=0 


(106) 


where 


Thus , 


x = cosa>. 


y = cosa>. 


V*> 

Tn(y) 


2 

cosco^m 


coso>2n 


A N 

K(x,y)  = Y x 
m=0 


m 


r n -1 

Y b(m,n)yn 
- n=0  1 


(107) 

(108) 
(109) 
(no) 

(111) 


where  m denotes  that  m is  fixed.  By  expanding  Eq  (111) 
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H(x,y)  = x°  J b(0,n)yn  + x1  £ b(I,n)yn 
n=0  n=0 


N 


+ . . . + xN  Y b(N,n)yn 
n=0 


(112) 


where  each  individual  summation  now  has  the  form  of  Eq  (99) • 
One  can  use  the  relation  defined  by  Eq  (105) *• 


H(x,y)  = x°  £ a(0,n)T  (y)  + . . . + xN  Y a(R,n)T  (y) 
n=0  n n=0  n 


N m N » 

= £ xm  * £ a(m,n)T  (y) 
m=n  n=0 


m=0 

N 


N 


= £ x • X a(m,n)cosw„n 

m=0  n=0 

n r N * -i 

= Y coscu^n  I y a(m,fi)xm 
n=0  ^ Lm=0  J 


(113) 

(114) 

(115) 

(116) 


where  the  term  in  brackets  in  Eq  (116)  is  also  of  the  same 
form  as  Eq  (99)-  Repeating  the  same  sequence  of  steps  yields 

N N 

H(x,w_ ) = £ cos<u?n  • y a(m,n)T  (x)  (11?) 

n=0  * m=0  m 

N N A 

= X cos<u0n  • y a(m,n)cosw,m  (118) 

n=0  z m=0  1 


The  notation  a(m,n)  has  no  mathematical  significance. 
It  is  used  only  to  emphasize  the  distinction  from  st(m,n). 
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(119) 


A N N A 

= X X a(m,n)coso>^m>  cosaj^n 

m=0  n=0 


which  is  the  desired  final  form. 

As  in  Chapter  II,  the  assumption  is  made  that  the 
impulse  response  is  centered  at  the  origin.  One  can  now 
compute  the  impulse  response  coefficients  {h(m,n)}  as 
follows : 

A A N NA_ 

= a(0,0)  + X a(m,  O)cosuj^m  + X a(0 ,n) cos^n 
m=l  n=l 


N N A 

+ X X a(m,n)ccsw,m  coso>pn  (120) 
m=l  n=l  ^ 


The  first  term  in  Eq  (120)  is  the  impulse  response  at  (0,0). 
The  second  and  third  terms  can  be  recognized  as  of  the 
same  form  as  Eq  (7)  of  Chapter  II.  The  fourth  term  must 
be  further  decomposed: 

N N A N N A 

X X a(m,n)cosw,m  cosw,n  = X cosG>,m  T a(m,n)costo~n  (121) 
m=l  n=l  x m=l  1 n=l  2 


N N * _ 

= X cosG>,m  X 2h(m,n)cosa>0n 
m=l  1 n=l  2 


11  * 

= 2cos&>,  X h(I ,n)cos<w_n  + ... 
x n=l  2 


” * 

+ 2cos(N&>,  ) X h(R  ,n)  cosw,n 


(122) 


(123) 
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• • • 


N A 

= 2cos&>,  2.  2h(I  ,n)cosa>2n  + 

n=l 

N a _ 

+ 2cos(Na.*,)  Y 2h(N,n)cosaj„n 

1 ri=l  * 

N N A 

= y y 4h(m,n)cosa>  m cosw„n 
m=l  n=l  1 2 


(124) 

(125) 


Thus,  knowing  {a(m,n)l , one  can  easily  compute  | h(m,n)} : 


h(0, 0)  = a( 0 , 0) , m=0 , n=0  (126) 
h(m,0)  = 1/2  a(m,0),  m/0,  n=0  (127) 
h(0,n)  = l/ 2 a(0,n),  m=0,  n/0  (128) 
h(m,n)  = l/4  a(m,n),  m/0,  n/0  (129) 

Comparison  of  Optimality 


It  is  sometimes  claimed  that  the  filters  designed  by 
the  McClellan  Transformation  are  truly  optimal.  McClellan 
made  this  claim  for  his  original  circularly  symmetric 
filters  and  later  extended  the  claim  to  include  "fan" 
filters  (Fig. 18).  McClellan  argued  that  along  any  slice 
of  the  2-D  filter,  the  approximation  problem  degenerates 
into  a 1-D  problem.  In  the  case  of  circular  symmetry  with 
t^  = .5,  t2  = = 0.5,  one  considers  the  slice  = 0. 

The  2-D  frequency  response  then  becomes  a 1-D  frequency 
response  ina»^« 

A 

F(&>^,0)  = -.5  + .5cos<u^  + .5  + .5  COSci>^  (130) 


(13D 


F(o>^)  = cosw^ 

Since  this  is  exactly  the  same  function  that  was  used  to 
generate  the  optimal  1-D  prototype,  the  2-D  filter  is 
claimed  to  also  be  optimal  (Ref  1:45-46).  Rabiner  and 
Gold  have  supported  this  contention  (Ref  7:477-478); 
however,  several  other  filtering  authorities  continue  to 
refer  to  the  McClellan  Transformation  as  sifff-optimal . 

Ideally,  empirical  tests  should  be  conducted  to 
measure  McClellan  Transformation  filters  against  known 
optimal  filters.  However,  the  design  algorithm  presented 
in  this  investigation  would  require  some  modifications. 

For  a fair  comparison,  both  of  the  2-D  transition  band 
edges  should  best  approximate  the  desired  contour  shape 
(as  is  done  in  true  optimal  methods). 

Table  II  records  data  obtained  on  the  design  of 
several  circularly  symmetric  optimal  lowpass  filters 
designed  by  Harris  and  Mercereau  (Ref  9:11.2).  Each  has  a 
passband  radius  of  0.4tr  and  a stopband  radius  of  0.6-rr. 

Using  the  original  McClellan  transformation,  the  0.4 
radius  corresponded  to  a 1-D  frequency  of  0.4000tt, 
while  the  0.6tt  radius  corresponded  to  a 1-D  frequency 
of  0.5760-rr  (recall  that  the  contours  become  non-circular 
as  a»-rr).  From  Table  II,  one  notes  that  the  McClellan 
error  is  somewhat  larger  than  the  optimal  error  for  each 
case.  If  the  McClellan  Transformation  is  used  with  1-D 
frequencies  of  0.4tt  and  0.6tt,  the  error  naturally  decreases. 
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Very  surprisingly,  it  actually  decreases  to  less  than  the 
optimal  error.  However,  this  second  case  is  not  a fair 
comparison:  it  would  be  better  to  design  the  optimal 
filters  with  respect  to  the  contours  resulting  from  this 
particular  McClellan  Transformation  and  then  compare  the 
errors.  Time  has  not  permitted  such  an  experiment  in 
this  investigation. 

The  data  available  suggests  that  the  McClellan 
Transformation  is  sub-optimal,  even  for  circularly 
symmetric  filters.  However,  the  absolute  difference  between 
the  optimal  error  and  the  McClellan  error  is  very  slight. 
Because  of  the  general  design  flexibility  and  significant 
savings  in  computation  time,  the  McClellan  Transformation 
compares  very  favorably  with  optimal  methods  for  designing 
2-D  FIR  filters. 


1 r 
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COMPARISON  OF  THE  ERROR  IN  OPTIMAL  AND  MCCLELLAN  FILTERS 


CDC  6600  Computer 


VI.  Design  Results 


This  chapter  presents  the  results  from  several  typical 
2-D  filter  design  problems  using  the  computer  program 
developed  in  this  investigation.  In  the  first  tests, 
simple,  familiar  filters  were  designed  to  validate  the 
design  program.  Gradually,  filters  with  more  novel 
specifications  were  designed  to  explore  the  limits  of  the 
method. 

Circularly  Symmetric  Filters 

To  test  the  design  program  against  a known  result, 
circularly  symmetric  filters  were  designed.  By  overriding 
the  default  constraint  and  inputting  Eqs  ( 62 ) — ( 66 ) , the 
original  McClellan  Transformation  was  duplicated.  The 
associated  contour  mapping  and  frequency  response  for  a 
21x21  point  filter  were  shown  in  Figs.  1 and  2,  Chapter  I. 
Fig.  10  illustrates  the  periodic  nature  of  2-D  digital 
filters . 

Circular  symmetry  was  also  obtained  using  the  default 
constraint.  Before  scaling,  the  range  of  the  mapping  was 
inadequate  (Fig.  11).  The  scaling  algorithm  was  used  to 
spread  out  the  contours , and  the  resulting  mapping  is  shown 
in  Fig.  12.  Comparing  Figs.  1 and  12,  it  is  seen  that  the 
contours  are  essentially  the  same  for  values  of  o>  up  to  about 
0.6tt.  (Extremely  close  comparison  reveals  that  the  contours 
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Fig.  11.  Circular  Contours  Using  the  Default  Constraint. 

w0  = .5tr,  tT  = [-.6377.  -.6377,  .6377,  .3622], 
5 = .00003 
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Fig.  12.  Circular  Contours  After  Scaling.  o>Q  = ,ky. 

tT  = [-.2840,  .5000,  .5000,  ,284o],  i = ,00i 

Default  constraint  with  scaling.  Prototype 
N = 21,  o>  = ,431xr,  « = .53tt,  Wt  = lil,  6 


of  Fig.  12  are  actually  more  circular.)  The  contours 
spread  out  differently  as  w — tt.  Essentially  identical 
frequency  responses  can  he  achieved,  as  shown  in  Fig.  13. 
However,  to  obtain  the  same  filter  specifications,  it  was 
necessary  to  use  1-D  prototypes  of  different  order.  Be- 
cause the  original  filter  was  designed  with  a lower  order 
prototype,  the  original  is  judged  superior. 

As  an  incidental  point,  note  that  the  information  gained 
from  the  frequency  response  drawings  is  also  available  from 
contour  mappings  of  the  extremal  frequencies  of  the  1-D 
prototype.  These  contours  precisely  define  the  ripple  in 
the  2-D  filter. 

Other  smooth  shapes  successfully  designed  included 
elliptical  shapes  (Figs.  14  and  15).  parabolic  shapes 
(Figs.  16  and  17).  and  hyperbolic  shapes  (Figs.  31  and  32 
of  Appendix  D) . 

Fan  Filters 

Fan  filters  are  used  in  seismic  analysis  to  filter 
unwanted  velocity  components  (Ref  26).  An  ideal  fan  filter 
is  shown  in  Fig.  18.  McClellan  developed  four  constraint 
equations  and  was  able  to  specify  t as 

tT  = [0.0,  0.5.  -0.5.  0.0]  (132) 

Fan  filters  were  designed  using  the  default  constraint, 
however  the  resulting  value  of  t 

tT  = [0.0003,  -0.4992,  0.5008,  -O.OOO3]  (I33) 
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13.  Comparison  of  Circularly  Symmetric  Filters. 

(a)  N = 19,  = ,5rr,  <ug  = .6tt,  Wt  = 1*1,  6=  .081, 

tT  = [-.5.  .5.  -5,  .5].  5 = .00000.  (b)  N = 21, 

w = .431tt,  o»s  = . 53tt  , Wt  = lil,  6=  .065, 

V - [-.284,  .5,  .5,  .284],  5 = .00003 
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W— RX I S (X  PI) 

Fig.  14.  Elliptical-Shaped  Filter  Contours.  Major  axis  = . 5tt 
Minor  Axis  = . 25tt.  N = 19,  o»  = . 2498rr,  w = .32rr, 

Pm  S 

Wt  = 1:2,  = .1 77,  62  = .08?,  tT  = [-.0482,  .7559, 

.0482,  .2440],  e = .00016. 
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Fig.  15.  Frequency  Response  of  Elliptical -Shaped  Filter. 
The  specifications  are  the  same  as  Fig.  14. 


Fig.  17.  Frequency  Response  of  Parabolic-Shaped  High  Pass  Filter. 
The  specifications  are  the  same  as  Fig.  17. 


reverses  the  sense  of  the  regions  defined  by  Fig.  18. 

One  could  remedy  this  by  using  a high  pass  1-D  prototype. 
Another  method  is  to  use  an  alternate  constraint,  0 — (tt.tt), 
and  then  scale  the  result.  This  produces  the  values  of 
Eq  (130) . The  resulting  contour  mapping  and  a frequency 
response  are  shown  in  Figs.  19  and  20. 

There  are  many  possible  variations  of  2-D  filters 
defined  by  straight  lines.  A diamond-shaped  filter  was 
specified  by  the  passband  contour 

a>2  = tt/2  - (134) 

The  resulting  mapping  and  a frequency  response  are  shown 
in  Figs.  21  and  22. 


Fig.  18.  Ideal  "Fan"  Filter 
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specifications  are  the  same  as  in  Fig.  19. 


Filters  with  More  Complex  Contours 


To  establish  the  limits  of  the  design  program,  several 
novel  filter  shapes  were  designed.  The  results  were 
generally  disappointing,  indicating  that  the  method 
cannot  deal  with  complex  contour  shapes  effectively. 

Oval-Shaped  Contours . The  desired  contour  for  an 
oval-shaped  filter  was  specified  as 


(a 


2 


tt/2 


01^  = 


'2  < Tf/2 

(135) 

»2>  tt/2 

(136) 

Eq  (135)  describes  a circle  of  radius  0.5tt  centered  on  the 
o>2-axis.  The  scaled  contour  mapping  is  shown  in  Fig.  23. 
As  tu  varies  either  side  of  0/ , the  contours  change  dramat- 
ically. In  this  particular  case,  an  assumption  made  in 
Chapter  IV  is  clearly  wrong:  the  immediate  neighbors  of 
this  contour  do  not  maintain  the  same  shape. 

Contour  Defined  by  a Cubic.  A contour  specificed  by 
the  cubic 


"2  = Vtt2  - rt/2p  + -rr/2 


(137) 


produced  a mapping  with  a very  large  average  error.  The 
contour  approximation  is  worst  in  the  vicinity  of  the 
* origin.  By  using  a properly  weighted  least-squares  criteria, 
the  fit  might  be  improved.  However,  the  primary  fault  is 
attributed  to  the  approximating  function,  Eq  (1).  The 
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Fig.  23*  Oval-Shaped  Filter  Contours.  o>Q  = . 502rr, 
tT  = [-.0174,  .0174,  .5121,  .4879  ], 

5 = .0047.  Default  constraint  with  scaling. 
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o 

o 


C 


Fig.  24.  Cubic  Contours.  a>Q  = .40774-rr,  tT  = [.2314,  -.2314, 
.6000,  .400oJ,  e = 0.2292.  Default  constraint 
with  scaling. 
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function  does  not  have  enough  terras  to  fit  complicated 
curves.  Adoption  of  higher  order  cosine  functions  in  a 
generalized  McClellan  Transformation,  Eq  (46),  would  solve 
this  problem. 

Non-monotonic  Contours.  In  Chapter  IV,  it  was  proven 
that  the  contour  mapping  must  be  monotonic.  All  empirical 
evidence  verifies  this;  a monotonic  contour  has  never  been 
observed.  This  does  not  imply  that  all  contours  for  a 
given  transformation  will  be  monotonic  in  the  same  direc- 
tion. However,  each  individual  contour  will  be  monotonic. 
All  attempts  to  specify  a non-monotonic  contour  led  to 
ill-defined  mappings  that  bore  no  resemblance  to  the 
desired  contour.  For  example,  the  semi-circle 


(Dr 


= ^-rr/2)2  - (a^  - tt/2)‘ 


(138) 


resulted  in  the  contour  mapping  shown  in  Fig.  25. 

Contours  not  Continuously  Differentiable . Several 
contours  that  were  not  continuously  differentiable  (non- 
smooth) were  used  to  test  the  design  program.  Immediately, 
problems  arose.  Fig.  26  shows  the  contours  resulting  from 
an  attempt  to  specify  the  right-angle  contour 


(Dr 


= Tf/2 , o»i  < rr/2 


w2  = °* 


> tt/2 


(139) 

(140) 


When  the  default  constraint  was  used,  the  mapping  was  not 
well-defined,  even  after  scaling.  By  imposing  the  addi- 
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25.  Non-Monotonic  Contour  Specification,  to  = .40342tt, 
tT  = [.1998,  -.1998,  .2577,  .7423],  5 = .1286 


Wl-flXIS  (X  PI) 


Fig.  27.  Right-Angle  Contours:  Alternate  Constraints. 

u>Q  = • 5tt , tT  =[-.2356,  .3399,  .6601,  .2356], 
e = .1465.  Constraint:  0— (0,0),  tr  — (tr , tt ) 


tional  constraint  tt  — (tt  , -rr ) , the  mapping  became  well- 
defined  but  the  average  error  was  large  (Fig.  27).  For 
non-smooth  increasing  contours,  the  best  constraints  to 
use  were  0 — ( 0 , tt ) and  it  — (tt,0).  These  constraints  were 
suggested  by  McClellan  (Ref  1:36-38),  but  they  do  not 
always  guarantee  that  the  range  of  the  mapping  will  be 
adequate.  However,  the  scaling  algorithm  can  still  be 
used  to  spread  out  the  contours. 


8^ 


VII.  Conclusions  and  Recommendations 


Conclusions 

The  primary  objective  of  this  investigation  was  to 
provide  a computer-aided  method  for  designing  2-D  digital 
filters.  Available  optimal  design  methods  were  rejected 
because  of  their  extreme  inefficiency.  The  McClellan 
Transformation  provided  an  effective  method  for  designing 
sub-optimal  filters  that  can  meet  a wide  variety  of 
specifications.  Design  time  for  the  filter  varies  in 
rough  proportion  with  the  square  of  the  filter  order 
(Fig.  28).  For  filters  exceeding  21x21  points,  the  compu- 
tational time  increases  rapidly.  Even  still,  the  design 
time  compares  very  favorably  with  other  available  methods. 
The  bulk  of  execution  time  is  spent  calculating  the  im- 
pulse response  coefficients;  if  these  are  not  required, 
the  design  time  is  shortened  considerably  (Fig.  28). 

The  coded  and  tested  FORTRAN  program  is  available  to 
Air  Force  Institute  of  Technology  users  on  the  CDC  6600 
computer.  The  structure  of  the  program  is  outlined  in 
Appendix  C,  and  the  use  of  the  program  is  explained  in 
Appendix  D.  A complete  design  example  is  included. 

The  limitations  of  the  method  must  be  recognized: 

1.  The  algorithm  maps  a specific  1-D  frequency  to  an 
approximation  of  a desired  contour.  The  "goodness"  of  the 
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contour  approximation  must  be  evaluated  in  each  design. 
Contours  with  high  curvature  or  comers  generally  lead 
to  unacceptable  mappings  due  to  large  rms  errors . 

2.  The  specified  contour  must  be  monotonic  in  the 
region  £o,Ttj  , limiting  the  flexibility. 

3.  For  some  contours,  the  default  constraints  may 
prove  inadequate.  The  suggestions  in  Chapter  VI  should 
be  considered. 


Fig.  28.  Approximate  Design  Times  for  2-D  Design 
Program 
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Recommendations  for  Future  Work 


Two-dimensional  digital  filtering  is  an  open  field 
for  research.  This  investigation  has  concentrated  on 
optimal  non-recursive  filters.  Recursive  filter  design  is 
worth  a full  investigation.  Little  has  been  said  regarding 
practical  realization  or  implementation  of  2-D  filters.  As 
an  aide  to  research  in  this  important  area,  a brief 
supplementary  bibliography  has  been  prepared.  Also  of 
importance  are  studies  of  the  practical  applications  of 
2-D  filters. 

There  are  several  areas  of  further  work  directly 
related  to  this  present  investigation.  Each  of  the 
following  suggestions  would  make  the  2-D  design  program 
more  useful  and  powerful » 

1.  It  would  be  advantageous  to  find  t such  that  two 
or  more  contours  are  approximated  simultaneously,  accord- 
ing to  an  error  criteria. 

2.  It  would  be  advantageous  to  adopt  the  general- 
ized McClellan  Transformation,  Eq  (46),  to  permit  better 
approximations  of  complicated  contours. 

3.  It  would  be  advantageous  to  calculate  the  im- 
pulse response  coefficients  by  a more  efficient  method. 

An  alternate  recursion  has  been  described  as  being  faster 
than  the  polynomial  expansion  used  in  this  investigation 
(Ref  26). 

ft  ! c 
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Appendix  A 


Algorithm  to  Determine  U for  Eg  (40) 


Given  the  three  Chebychev  polynomial  identities 


Tq(x)  = 1 

(141) 

Tx(x)  = x 

(142) 

Wx)  = 2xTn(x)  - VjIx) 

(143) 

one  can  write  an  expression  for  any  Tm(x).  Each  of  these 
equations  could,  in  turn,  be  used  to  solve  for  the  powers 
x°,  x^,...,xn.  Table  III  illustrates  the  first  few  terms 
of  such  a listing. 

A recursive  algorithm  for  determining  the  coefficients 

O 1 

of  Tq(x),  T1(x) , . . . ,Tn(x)  in  each  expression  for  x , x , ...x 
has  been  developed.  The  algorithm  was  deduced  by  writing 
the  equations  for  x1  in  matrix  forms 


1 

0 

1/2 

0 

3/8 


0 0 

1 0 

0 1/2 

3/4  0 

0 4/8 


0 0 

0 0 

0 0 

1/4  0 

0 1/8 


Tq(x) 

T-^x) 

T2(x) 

T3(x) 

V*) 


(144) 


Tn(x) 
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TABLE  III 


The  Chebychev  Polynomials  Solved  for  x 


T (x)  = 1 = x° 

0 

Tx(x)  = x 

x°  = T0 
x1  = T2 

T2(x)  = 2x2  - 1 

x2  = 1/2 

*0  + *2  J 

T^(x)  = 4x3  - 3x 

x3  = 1/4 

3T1  + *3. 

\ 

T4(x)  = 8x4  - 8x2  + 1 

• 

• 

• 

x4  = 1/8 

• 

• 

• 

3T  + 4 T. 

• 0 i 

2 + **] 

If  the  rows  of  Eq  (144)  are  partitioned, 


x°  = 

[1 

0 

0 

0 

0 • • • 

1-4 

T 

(145) 

x1  = 

[° 

1 

0 

0 

0 • • • 

] = 2 !* 

T 

(146) 

x2  = 

1/2  [1 

0 

1 

0 

0 • • § 

1=2* 

T 

(147) 

X3  = 

1/4  [ 0 

3 

0 

1 

0 . . . 

1=2* 

T 

(148) 

X4  = 

1/8  [3 

0 

4 

0 

1 

1=2? 

T 

(149) 

it  is  noted  that  for  m > 1 


(150) 


where 


Um, 0 " Um-1, 1 


(15D 
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(152) 


U , = 2U  , + um  , 

m, 1 m-l, 0 m-l,2 


Um,k  Um-l,k-l  + Um-l,k+l 


U . = 1 

m,k 


U > = 0 
m,  k 


k<m 

(153) 

k=m 

(154) 

k>m 

(155) 

\ 

Thus,  a recursive  formula  for  each  expression  x°,  x1,...^11 
is  available : 

xm  = l/2m_1  U*  T (156) 


c 
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Appendix  B 


The  Method  of  Least-Squares 
Curve  Fitting  with  Constraints 

Case  I:  No  constraints  on  t. 

A known  (or  measured)  function  y = f(x)  is  given.  It  is 
desired  to  approximate  y with  the  function 


Pm(x)  = tlSl(x)  + t2g2(x)  + ...  + tmgm(x) 


(157) 


where  each  g^(x)  is  linearly  independent.  Defining 


6k  “ pm^xk^  ~ yk’  k " 1,2 N 


(158) 


it  is  desired  to  minimize 


A(i)  = Jo  6 J = Jo  IV**’  ' yJ: 


(159) 


where  the  coefficients  t are  unknown  and  unconstrained.  Let 


' V = *kj  = J0  gi(xi)gk<xi> 


(160) 


“ll  “l2 


- = “21  “22 


(161) 


X yigi(xi} 

1=0  j 


(162) 


?T  = [»! 


(163) 
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It  can  be  shown  that  (Ref  27  022-323) 

« t = £ (16*0 

or 

t = of1  £ (165) 

is  the  best  least-squares  fit  of  the  approximating  function 
p (x)  to  the  data  set  y,  provided  that  a is  non-singular. 

Case  II : One  constraint  on  t. 

The  constraint  on  t has  the  form 

h(t)  = c-,t,  + ...  + c t - c = 0 (166) 

- i±  m m o 

A necessary  condition  for  an  extremum  is  that  the  total 
differential*  vanish: 

pt  dtl  + + pt  dtm  = 0 (16?) 

1 m 

The  total  differential  of  the  constraint  is 

h.  dt,  + . . . + h.  dt  =0  (168) 

t -i  1 t m 

1 m 

If  Eq  (168)  is  multiplied  by  \ , a Lagrangian  multiplier, 
and  added  to  Eq  (167)  (Ref  28:169-170) 


i ( p^  is  shorthand  notation  for  3p/at1. 
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0 


(169) 


■4 


(170) 


P+  + =0  (171) 

m m 

h(t)  = 0 (172) 

this  gives  (m  + 1)  equations  in  the  m + 1 unknowns:  t,  X.. 

As  in  Case  I,  Eq  (159)  for  A(t)  can  be  written.  If  one 
multiplies  the  term  in  brackets  and  performs  a partial 
differentiation  with  respect  to  t, 

N 2 

At1  = i^0^2tlgl  + 2t2glg2  + ~ 2giyi^  (173) 

t 

■ 

Atm  = ilo(2tlglgm  + + 2tmgm2  “ 2gmyi)  (1?4) 

Because  Eqs  (170)  - (172)  were  determined  for  any  p^,  one  can 
substitute  A^  for  p^.  Letting 


N 


results  in 


By  noting  that 


and  that 


■«  - 1/2  s« 


= 1/2  I, 


2h(t ) = 0 


(177) 

(178) 

(179) 

(180) 

(181) 

(182) 

(183) 


is  an  equivalent  set  of  constraints 


h 

• 

• 

• 

il 

o 

o 

‘11 


r21 


aml 

Ci 


*lm 

*2111 


'mm 


m 


it  is  easily  shown  that 


(184) 
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* /*x-l  * 

t = («)  1 £ 


so  long  as  a is  invertable. 


(185) 


Case  Ills  Multiple  constraints  on  t. 

If  there  are  n<m  constraints  on  t,  there  will  be  n 


equations  of  the  form 


h.  (t ) — c - t -|  + •••  + c.  t — c . — 0 

l n i lm  m io 


(186) 


For  each  constraint,  a Lagrangian  multiplier  is  introduced, 
resulting  in  the  following  set  of  equations: 


Atx  + Xlhtx  + • • • + Xnht^  - 0 


At  + ^ht  + . . . + Xnht  - 0 
mm  m 


h1(t)  = 0 


hn(t)  = 0 


(187) 


(188) 


(189) 


(190) 


This  gives  (m  + n)  equations  in  the  (m  + n)  unknowns:  t, 

Xl’  X2’ * * * ’ Xn* 

Extending  the  argument  of  Case  II,  it  can  be  shown  that 


99 


Appendix  C 


Comments  on  Programming 

Program  Structure 

In  structuring  this  program,  the  primary  goals  were  that 
the  program  would  be  understandable  and  modifiable  by  subse- 
quent users.  To  permit  this,  the  principles  of  modem  soft- 
ware engineering  have  been  followed  as  outlined  in  the  text 
by  Yourdan  and  Constantine  (Ref  29). 

At  the  start  of  this  investigation,  the  program  re- 
quirements were  defined: 

1.  The  program  would  be  interactive.  This  would 
provide  maximum  user  flexibility.  Additionally,  it  would 
allow  the  program  to  be  semi-educational,  potentially  useful 
in  a course  on  digital  filtering. 

2.  The  program  would  design  two-dimensional  filters 
with  arbitrary  transition  band  contours.  (It  was  later 
discovered  that  there  is  an  admittable  class  of  contours.) 

3.  The  minimum  output  would  include  the  array  of  2-D 
impulse  response  coefficients  and  the  frequency  response  of 
the  filter.  There  would  be  an  optional  provision  for  a 
computer  plot  of  the  frequency  response. 

4.  The  program  would  be  simple  to  use,  not  requiring 
a sophisticated  knowledge  of  SCOPE  procedures  (Ref  32). 

A high  level  structure  chart  of  a first-cut  design  is 
shown  in  Fig.  29  Each  of  the  modules  represents  an 
independent  process . 


Fig.  29.  First-Cut  Structure  Chart  for  Design  Program 


Considerations  for  Interactive  Use 

Because  the  program  was  required  to  be  interactive, 
memory  management  had  to  be  considered.  Using  Control  Data 
Corporation  INTERCOM,  local  procedures  establish  a maximum 
field  length  of  60K  octal  memory.  This  prevents  a large 
program,  such  as  this  one,  from  being  loaded.  There  are 
two  methods  available  to  circumvent  this  restriction; 
segmentation  and  overlays  (Ref  30:1-12).  Although  segmen- 
tation is  a more  elegant  method,  overlays  were  selected  to 
permit  an  eventual  interface  with  a large  package  of  computer- 
aided  design  programs  concurrently  being  assembled  (Ref  31)* 
Overlays  may  extend  three  levels  deep.  To  permit  an  inter- 
face with  another  program,  one  overlay  level  must  be  left 
for  communication,  leaving  only  two  levels  to  work  with. 

The  main  level  overlay  has  been  used  to  serve  as  a 
controller  for  the  design  program.  It  corresponds  to  the 
top  box  in  the  structure  chart  of  Fig.  29.  The  initial  plan 
was  to  assign  each  of  the  other  modules  to  a separate  pri- 
mary level  overlay.  Unfortunately,  with  this  program 
division,  the  60K  restriction  could  not  be  met.  Thus, 
several  of  the  overlays  have  been  artificially  divided  into 
two  distinct  primary  overlays , with  a great  deal  of  addi- 
tional control  passing  required.  This  compromise  hurt  the 
program  structure  by  increasing  the  dependence  of  the 
modules,  but  it  was  necessary  to  permit  interactive  use. 

o 
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Avoiding  SCOPE  Commands 

Another  troublesome  aspect  of  the  program  involves  the 
"Get  Functional  Contour"  module  of  Fig.  29.  The  program 
was  required  to  design  filters  with  arbitrary  transition 
band  contours  (within  the  admittable  class  of  contours), 
thus  it  was  impossible  to  develop  a list  of  available 
functions  for  the  user.  The  only  realistic  approach  was  to 
allow  the  user  to  create  his  own  subroutine  to  describe 
, the  contour  shape.  This  subroutine  is  trivial  to  write; 

however,  integrating  such  a subroutine  with  the  binary  code 
of  the  design  program  is  difficult.  The  subroutine  must  be 
compiled  and  libraried,  a tedious  task  involving  a good 
knowledge  of  SCOPE  commands  (Ref  32). 

Fortunately,  a method  was  found  that  relieves  the  user 
of  any  undue  burden  while  still  permitting  him  to  supply 
the  required  subroutine.  A "cataloged  procedure”,  written 
using  University  of  Washington  Control  Language  (Ref  33), 
performs  all  the  functions  associated  with  compiling  and 
librarying  the  subroutine,  attaching  and  executing  the 
design  program  (MCLTRN),  and  routing  a hardcopy  of  the  out- 
put listing  and  plots.  The  user  creates  his  subroutine, 
saves  it  as  a local  file,  and  then  issues  a set  of  simple 
commands.  The  nature  of  the  subroutine  and  the  required 
procedure  commands  are  outlined  in  detail  in  Appendix  D. 

A listing  of  the  cataloged  procedure  (PROFIL)  is  attached 
to  this  appendix. 

C 
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Program  Flowchart 

This  program  was  coded  with  the  concepts  of  structured 
programming  in  mind.  Unconditional  branching  has  been 
avoided  whenever  a reasonable  alternative  existed.  Sub- 
routine calls  are  made  very  freely;  however,  most  subroutines 
are  short  and  functional.  As  a rule  of  thumb,  modules  con- 
tain between  10  and  50  lines  of  code.  For  these  reasons, 
detailed  flowcharts  for  each  module  are  unnecessary.  Fig.  30 
shows  a top-level  flowchart  for  the  design  program.  The 
flowchart  has  been  deliberately  held  to  one  page , showing 
only  the  detail  required  for  a quick  overall  view  of  the 
program. 

Program  Source  Code 

The  design  program  source  code,  PROGRAM  MCLTRN;  (Ref  35) 
involves  over  1600  statements,  so  it  has  not  been  included 
in  this  report.  Copies  of  the  source  listing  are  available 
from  the  author  or  from  Prof.  Gary  Lamont,  AFIT/EN. 

Program  PROFIL 

The  cataloged  procedure  PROFIL  controls  Program  MCLTRN . 

A listing  follows: 
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Program  PROFIL 


COMPILE (NAME=X) 

RETURN, GOGO. 

REWIND, NAME. 

FTN , I=NAME , L='0,  B=GOGO  . 

EDITLIB , I=DA , L=LIST . 

RETURN, LIST, DA. 

LIBRARY, MYLIB. 

/DATA  DA 

LIBRARY (MYLIB, NEW) 

ADD (*, GOGO) 

FINISH. 

ENDRUN . 

RECOMPILE (NAME=X ) 

RETURN, GOGO. 

REWIND, NAME. 

FTN , I=NAME , L=0 , B=GOGO . 

EDITLIB , I=DB , L=LIS  T . 

RETURN, LIS T,DB. 

LIBRARY, MYLIB. 

/DATA  DB 

LIBRARY (MYLIB  ,OLD) 

DELETE (*) 

ADD (*, GOGO) 

FINISH. 

ENDRUN . 

DESIGN 

IF ( -FILE , AFITSUB ) 

ATTACH, AFITSUBROUTINES , ID=AFIT . 
LIBRARY , AFITSUB , MYLIB . 

IF (-FILE, MX) 

ATTACH , MX , THES IS , C Y=3/2 , ID=RBB . 
REWIND, MX. 

MX. 

RETURN, MCL. 

ROUTE ( TERM=BB , USER=K J/ ) 

REWIND, RESULT. 

COPYSBF , RESULT , RRX . 

ROUTE , RRX , DC=PR , TID=TERM , FID=USER . 
RETURN , RESULT . 

IF (FILE, PLOT) 

ROUTE , PLOT , DC=PT , TID=TERM, FID=USER . 


FILES. 

PLFILE 

REWIND, TAPE2. 

REQUEST , TEMP3D , *PF . 

COPY, TAPE2 , TEMP3D. 

CATALOG , TEMP3D , DATA3DPLOT . 

RETURN , TEMP 3D . 

IF ( -FILE , PLOT3D ) 

ATTACH , PLOT3D , THES IS , CY =4// , ID=RBB. 
REWIND, PLOT3D. 

BATCH , PLOT3D , INPUT , HERE . 

FILES. 

3DPLOT 

IF (-FILE, PLFILE) 

ATTACH , PLFILE , YOURFILE . 

IF (-FILE, DISSPLA) 

ATTACH, DISSPLA , ID=X654321 . 

LIBRARY, DISSPLA. 

ONLINE. 


NOTE  1 Separate  procedures  with  end-of-record. 


( 
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Appendix  D 


User's  Guide  to  Program  MCLTRN 


MCLTRN 

15  December  1977 

Identification : 

MCLTRN  - McClellan  Transformation  Design  for  Two- 
Dimensional  (2-D)  Digital  Filters 

FORTRAN  Extended  (FTN)  Program 

Air  Force  Institute  of  Technology 

Wright-Patterson  Air  Force  Base,  Ohio 

Capt.  R.  B.  Brown,  GE-77D 


Purpose « 

MCLTRN  will  design  2-D  finite  impulse  response  (FIR) 
zero-phase  filters  of  the  form 

N N 

G(u>i , cog ) = Z S a(m,n)cosw1m  cosw2n  (45) 

. m=0  n=0 


The  filter  will  be  a very  close  approximation  to  an  optimal 
filter.  Program  output  includes  a 2-D  frequency  response 
(optionally  plotted),  the  2-D  impulse  response  coefficients, 
and  contour  maps.  MCLTRN  is  intended  for  interactive  use. 


Control! 

ATTACH,  PROFIL,  ID  = T770314. 

A cataloged  procedure  (Ref  33)  PROFIL  is  used  to 
control  MCLTRN.  There  are  six  commands: 


BEGIN,  COMPILE. 
BEGIN,  RECOMPILE. 
BEGIN,  DESIGN. 


BEGIN,  ROUTE. 
BEGIN,  PLFILE . 
BEGIN,  3DPLOT. 


Programming  Information: 


1.  MCLTRN  requires  a user  supplied  FORTRAN  SUBROUTINE 
FUNCT(I,  N,  X,  Y),  where  X,  Y are  the  spatial  axes,  N is  the 
number  of  sample  points  (supplied  by  MCLTRN),  and  I is  the 
current  sample  index  (supplied  by  MCLTRN).  The  subroutine 
defines  a smooth  monotonic  contour  in  the  2-D  region 
[O.lTx  [0,1]  . MCLTRN  will  scale  this  contour  to  the  region 
LO , -rr  J x [0 , tt]  . 

e.  g.  To  describe  a quarter  circle,  centered  at  (0,0)  with 
radius  0 . 5» 

SUBROUTINE  FUNCT(I,  N,  X,  Y) 

C DEFINE  VALID  X RANGE  ON  [0,l]  AND  INCREMENT  X 
XMAX  = 0.5 

X = XMAX  * (I  - 1.0)  / N 
C CALCULATE  CORRESPONDING  Y VALUE 
Y = SQRT  (XMAX**2  - X**2) 

RETURN 

END 


2.  After  creating  FUNCT  (I,  N,  X,  Y),  the  edit  file 
must  be  saved,  without  sequencing,  as  local  file  X: 

SAVE,  X,  NOSEQ,  OVER. 

To  automatically  compile  and  library  the  subroutine,  command: 
BEGIN,  COMPILE. 

If  there  are  FORTRAN  errors,  or  if  it  is  necessary  to  compile 
another  subroutine,  use  the  alternate  command: 

BEGIN,  RECOMPILE. 

(If  error  messages  related  to  your  user  library  MYLIB  appear, 
check  that  you  saved  file  X and  then  attempt  to  use  the 
recompile  command.) 

3.  To  execute  MCLTRN,  command: 

BEGIN,  DESIGN. 

4.  To  obtain  a hardcopy  of  the  output  listing  (RESULT) 

or  contour  plots  (PLOT),  command:  1 j 

( BEGIN,  ROUTE  (,  ,XX,  YYYYY). 
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where  the  grouping  in  parenthesis  is  optional: 

XX  = TID  code  (optional  default  = BB) 

YYYYY  = FID  code  (optional  default  = USER) 

Note:  Plots  will  always  have  the  user  identification 
code  as  the  banner. 

5.  To  create  a DISSPLA  file  for  the  3-D  drawing  of  the 
frequency  response,  command: 

BEGIN,  PLFILE . 

Because  DISSPLA  (Ref  34)  requires  over  100K  memory,  the 
job  will  be  batched  to  the  input  queue.  After  allowing  time 
for  execution,  check  your  output  file  JJJucxx  for  errors 
(uc  is  your  user  identification  code): 

BATCH,  JJJucxx,  LOCAL. 

EDIT,  JJJucxx,  SEQ. 

Assuming  there  were  no  FORTRAN  errors*,  you  may  create  (or 
preview)  and  dispose  a Calcomp  plot  file  by  commanding: 

BEGIN,  3DPL0T. 

DRAW  = 1-END$ 

If  you  are  at  a graphics  terminal,  command: 

TEK4010 

TEK4014  as  appropriate 

etc . 

To  dispose  your  plot,  command: 

ROUTE,  PLOT,  TID  = XX,  FID  = YYYYY,  DC  = PT. 


* 

In  particular,  FORTRAN  error  103:  No  permanent  file 
space  for  plot  file  — job  killed. 


Methods 


MCLTRN  maps  each  frequency  of  a 1-D  bandpass  filter  to  a 
contour  in  the  spatial  plane  by  the  McClellan  Transformation: 

cos o>  = t1  + t2  cos(x)  + Xj  cos(y)  + t^  cos(x)  cos  (y)  (192) 

In  particular,  the  frequency  coQ  (usually  a transition  band 
boundary)  is  mapped  to  the  contour  defined  by  SUBROUTINE 
FUNCT  (I,  N,  X,  Y).  MCLTRN  optimizes  t by  the  least-squares 
criteria.  For  a non-trivial  mapping,  constraints  on  t are 
necessary;  however,  a default  constraint  that  ft  maps  to 
OM)  is  available.  An  optional  contour  scaling  algorithm 
returning  a new  1-D  frequency  coQ,  is  recommended  to 
guarantee  a well-behaved  frequency  response.  The  user  may 
examine  the  locus  of  contours  generated  by  t as  a function 
of  eo,  in  both  tabular  or  plotted  form. 

Once  t is  established,  a 1-D  prototype  optimal  FIR 
zero-phase  filter  is  designed.  This  filter  is  normally  a 
low-pass  or  high-pass  filter  with  a stopband  or  passband 
frequency  of  coq . At  this  writing,  the  filter  must  be  of 
odd  order  not  exceeding  21.  The  error  magnitude  in  the 
1-D  prototype  will  be  passed  to  the  2-D  filter,  where  each 
1-D  extremal  frequency  maps  to  a 2-D  contour.  A plot  file 
of  this  mapping  is  created  automatically  by  MCLTRN. 

Finally,  MCLTRN  calculates  the  2-D  frequency  response. 
Optionally,  the  2-D  impulse  response  coefficients  will  be 
calculated. 


- 


Ill 


Special  Conventions : 


1.  The  2-D  filter  designed  using  the  McClellan 
Transformation  is  symmetric  about  both  the  x and  y axes. 

Thus,  only  a first  quadrant  design  is  considered. 

2.  All  frequencies  are  normalized  to  the  interval 

[0,Tf]  and  must  be  entered  as  multiples  of  tt  to  l.fi). 

Design  Example* 

The  following  example  illustrates  the  use  of  PROFIL 
and  MCLTRN.  The  desired  contour  is  defined  by  the 
rectangular  hyperbola 

xy  = a2/2  (193) 

where  a = 0.30tt.  The  resulting  contour  mapping  is  shown 
in  Fig.  31-  The  frequency  response  is  shown  in  Fig.  32. 

Note  that  user  entries  are  shown  in  lower-case  letters, 
programming  prompts  in  upper-case  letters. 
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19  X 19  POINT  FILTER 
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(after  allowing  time  for  execution) 
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Fig.  32.  Frequency  Response  of  Hyperbolic-Shaped  Filter. 
The  specifications  are  the  same  as  Fig.  31. 


Appendix  E 


Use  of  Advanced  DISSPLA 


DISSPLA  is  a software  graphics  package  developed  by- 
integrated  Software  Systems  Corporation.  The  DISSPLA 
Reference  Manual  (Ref  3M  is  augmented  locally  by  prelimi- 
nary and  post-processing  instructions  available  from  the 
ASD  Computer  Center. 

One  of  the  features  of  Advanced  DISSPLA  is  a capability 
for  projecting  three-dimensional  (3-D)  surfaces  and  lines 
(Ref  3^*Adv-C).  3-D  surfaces  are  plotted  with  "hidden 

lines"  automatically  removed.  Any  viewpoint  not  actually 
within  the  surface  being  drawn  is  allowed. 

The  basic  plotting  subroutine  for  3-D  surfaces  is  CALL 
SURFUN(ZFUN,  IXPTS , XDELTA , IYPTS,  YDELTA , WORK),  where 

ZFUN  is  the  function  name  for  a two-variable 
function  to  be  evaluated,  producing  a 
magnitude  z . 

IXPTS  is  the  number  of  sample  points  to  be 
evaluated  in  each  XDELTA  interval. 

IYPTS  is  analogous  to  IXPTS , but  over  the  YDELTA 
interval . 

XDELTA  is  the  interval  between  lines  parallel 

to  the  x-axis.  The  smaller  this  interval, 
the  greater  the  detail. 

YDELTA  is  the  interval  between  lines  parallel  to 
the  y-axis. 

WORK  is  an  internal  array,  dimensioned  at  least 
equal  to  the  number  of  perimeter  sample 
points  plus  4. 

c 
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The  limits  on  x,  y,  and  z are  specified  by  calls  to  GRAF3D 
and  AXES3D.  For  information  on  these  and  other  calls,  the 
DISSPLA  Advanced  Reference  Manual,  Part  C,  should  be  re- 
viewed. 

A source  listing  of  the  FORTRAN  program  used  to  generate 
the  3-D  frequency  response  drawings  of  this  investigation 
has  been  included  as  an  illustration  of  3-D  DISSPLA. 

Function  G evaluates  Eq  (44)  over  a grid  of  points  in  x-y. 
Tape  2 contains  t,  and  the  coefficients  (b(n)}. 


J J J , T20  , CP 120000.  T770314 . NCLTRN  USER 
PTN(R»0  * 

attach, tapes, cata3cplot . 

REt'lHD,TAFE£. 

ATT *CH, DISSPLA, I D-X2S4321. 

USRhR',.  DISSPLA. 

REOUEST , PLF I LE . »pF . 

LOO. 

CATALOG, PLFILE.VOURFrLE. 

PURGE, TAPE2. 

•COR 

.PROGRAM  PLT3D* INPUT, OUTPUT. TAP£2,PLFILE*0> 

DIMENSION  U0R):«6Se>.l<ll>.Ti4) 

common  /one-  b.np.t 

EXTERNAL  G 
READ* 2, 11 > NP.T 

11  FORMAT ( 1 3 , 4E 1 3 . 9 > 

READ'  2, 12MB<f1',N*l,NP) 

12  F0RNAT15E18.  »•  5E12.9--E18.9> 

PRINT  12, T 

print  i2,(B‘.M).h*i,nP) 

CALL  COMPRS 
CALL  BGNPL'.l) 

CALL  TITL30UH  .-1.7 .9.6.9) 

CALL  UUhBS' 12. ,10. .12. ) 

CALL  AXES3D'0.0.9, 0,9, 0,2. 0.2. ®,1.5) 

CALL  GRAF3D(-1. 0.0. 1,1.0. -1.0, 0.1. 1.0, -.25.0. 1,1. 25) 
CALL  SUPFUNIG.  3,. 125.  3,. 125.  UORK ) 

CALL  ENBPL'l) 

CALL  DOHEPL 

STOP 

END 


FUNCTION  G(U1,U2) 

COUPON  /ONE / B.NP.T 
DIMENSION  B(ll).T'4> 

PI *3. 1415926536 
UU1*U1«PZ 
UU2*U2lPI 
X-COS(UUl) 

V*C0S(LU2 ) 

2-T(4)*X»v 

X-T(2)*X 

V-T<3)»V 

Ki.m 

30  M-l.NP 

30  0*0  ♦ KN)*( CT(l >4X4Y*Z1M<P-1 )) 

RETURN 
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